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Abstract 


A  model  for  the  formation  and  evolution  of  three-dimensional  sedimentary  structures 
such  as  longshore  sand  ridges,  on  the  continental  shelf  in  water  deeper  than  that  of  the 
shoaling  region,  is  proposed.  The  model  is  based  on  the  interaction  between  surficial  or 
internal  weakly  nonlinear  shallow  w'ater  waves  having  weak  span-wise  Spatial  dependence 
and  the  bottom  topography. 

While  these  ridges  are  not  the  result  of  a  single  formative  agent,  it  is  argued  that  the 
mechanism  proposed  in  this  study  does  contribute  significantly  to  their  generation  and 
evolution.  Testing  the  hypothesis,  however,  is  as  difficult  as  formulating  it.  Comparisons 
of  this  model  with  oceanographic  data  must  wait  for  sufficient  data  to  become  available. 
In  conjunction  with  developing  the  sand  ridge  model,  this  study  proposes  a  new  math¬ 
ematical  equation,  properties  of  which  are  explored  here  in  some  detail.  This  equation 
potentially  applies  to  other  physical  processes  and  raises  questions  which  are  themselves 
good  avenues  for  further  research. 

The  numerical  implementation  of  the  model  combines  fixed  point  methods  with  finite 
difference  techniques,  resulting  in  a  scheme  which  is  found  to  be  superior  to  conventional 
finite  difference  techniques  in  economy  of  computational  resources  and  speed.  Details  of 
the  scheme’s  inner  workings  and  its  performance  are  included. 
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Chapter  1 

Introduction 


1.1  Statement  of  the  Problem 

The  dynamics  of  sand  ridges,  which  are  a  common  feature  of  the  Continental  Shelf, 
are  poorly  understood.  Sand  ridges  are  underwater  bar-like  features  composed  of  loose 
granular  sediment.  Hundreds  of  meters  long  and  up  to  a  few  meters  high  sand  ridges  are 
usually  found  in  groups,  arranged  in  more  or  less  parallel  rows  separated  from  each  other 
by  hundreds  of  meters.  They  may  be  loosely  classified  as  either  tidal  ridges  or  longshore 
sand  ridges.  Tidal  ridges  are  oriented  parallel  to  the  prevailing  direction  of  the  local  ocean 
currents,  whereas  longshore  sand  ridges  are  oriented  normal  to  the  direction  in  which 
the  overlying  water  waves  propagate.  In  this  study  I  propose  a  possible  mechanism  for 
the  formation  and  evolution  of  longshore  sand  ridges. 

The  model  presented  here  follows  from  work  initiated  by  Boczar-Karakie  .vicz  and 
Bona,  which  dates  to  1986.  In  [1]  they  conjecture  that  longshore  sand  ridges  are  the 
result  of  energetic  interactions  between  shallow  water  waves  and  the  underlying  bottom 
topography,  and  propose  a  simple  model,  which  in  [2]  was  shown  to  be  in  qualitative 
agreement  with  oceanic  data.  While  the  present  study  owes  much  to  the  previous  work, 
it  improves  upon  it  considerably  and  in  several  ways.  In  addition  to  extending  the  two- 
dimensional  model  to  three  dimensions,  this  work  contributes  to  an  understanding  of 
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the  general  behavior  and  mathematical  structure  of  both  the  two-  and  three-dimensional 
model.  The  model  is  free  of  adjustable  parameters  and.  at  this  stage,  intentionally  crude. 
Our  motivation  was  to  present  the  simplest  possible  formulation  in  order  to  effectively 
study  and  test  the  hypothesis  unhindered  by  physically  negligible  effects.  We  do  not 
claim  that  sand  ridge  formation  is  the  result  of  a  single  event  or  agent,  nor  do  we  claim 
that  this  model  rules  out  all  other  explanations  for  the  phenomenon.  Rather,  we  describe 
a  likely  mechanism  for  the  formation  and  evolution  of  these  structures,  a  mechanism  we 
believe  must  play  a  significant  role. 

A  great  deal  of  work  on  the  problems  of  sedimentation  has  been  done;  however,  par¬ 
ticularly  since  the  middle  of  this  century,  most  of  the  work  has  been  directed  towards 
understanding  smaller-scale  aspects  of  sediment  motion,  rather  than  the  formation  and 
evolution  of  sedimentary  structures.  For  a  comprehensive  review  of  the  present  level  of 
our  understanding  of  sedimentation,  the  reader  is  referred  to  Sleath  [3],  While  consid¬ 
erable  progress  has  been  made,  our  current  understanding  of  sediment  dynamics  and, 
especially,  of  sedimentary  structure  formation  is  far  from  complete.  The  emphasis  in  this 
dissertation  will  be  on  the  fluid  mechanical  aspects  of  sedimentation.  We  believe  that  a 
great  deal  of  progress  in  understanding  sediment  movement  in  a  fluid  environment  can 
be  achieved  by  determining  first  what  sort  of  patterns  the  fluid  is  able  to  generate. 

The  plan  of  the  dissertation  is  as  follows.  In  this  introduction,  we  discuss  the  rele¬ 
vance  of  the  study,  describe  the  morphology  of  oceanic  sedimentary  structures,  comment, 
on  observational  and  laboratory  work,  and  review  the  various  sedimentation  and  sand¬ 
bar  models.  In  Chapter  2,  we  consider  the  main  hydrodym  nical  issues  for  the  cases  of 
both  surface  and  internal  waves.  Chapter  3  deals  with  the  boundary-layer  problem  and 
with  the  development  of  a  mass  transport  equation  driven  by  the  nonlinear  wave  motion. 
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Chapter  4  presents  analytical  results  pertinent  to  the  new  equations  resulting  from  the 
hydrodynamic  problem  discussed  in  Chapter  2.  Chapter  5  presents  the  numerical  solu 
tion  of  the  full  system,  along  with  an  analysis  and  evaluation  of  the  numerical  scheme. 
Numerical  examples  are  presented  and  qualitatively  discussed  in  Chapter  6'  Chapter  7 
lists  conclusions  and  open  questions  worthy  of  future  pursuit.  Two  appendices  provide 
details  on  the  higher  order  theory. 

1.2  Relevance  of  This  Study 

For  many,  celebrating  the  beauty  and  mystery  of  nature  is  sufficient  reason  for  studying 
the  patterns  and  structures  by  which  nature  organizes  and  evolves.  Nevertheless,  there 
are  also  very  practical  reasons  for  research  into  sedimentary  structures,  some  of  which 
are  listed  below: 

•  The  study  and  control  of  coastal  erosion  is  of  major  economic,  political,  and  eco¬ 
logical  importance  to  communities  that  neighbor  oceans  and  major  lakes. 

•  Most  features  of  the  ocean  bottom  evolve  in  geological  time  scales;  sedimentary 
structures,  however,  change  comparatively  quickly.  A  model  of  sedimentary  struc¬ 
ture  evolution  and  movement  will  help  us  understand  how  these  quickly  evolving 
features  will  modify  the  bottom  topography  over  time  scales  relevant  for  such  things 
as  navigation. 

•  Understanding  the  movement  of  these  structures  may  help  biologists  discern  how- 
nutrients  and  organic  materials  migrate  along  the  ocean  bottom,  information  es¬ 
sential  to  understanding  the  dynamics  of  the  marine  habitat. 

•  Similarly,  such  knowledge  may  shed  light  on  the  movement  and  eventual  fate  of 
man-made  pollutants  and  debris. 
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•  Alternatively,  the  sedimentary  structures  themselves  may  have  economic  and  social 
importance.  The  best  surfing  beaches  have  naturally  occurring  sandbars  strategi¬ 
cally  located  to  concentrate  the  action  of  water  waves  in  some  areas  and  destroy  it 
in  others.  Predictability  of  these  so-called  "hot  spots"  is  essential  to  the  welfare  of 
the  surfing  community. 

•  Sand  ridges  are  part  of  hydrocarbon  reservoirs  in  ancient  strata.  Predicting  their 
properties  and  evolution  would  be  useful  in  petroleum  exploration. 

•  The  storm-wave  devastation  of  coastal  communities  and  offshore  structures  could, 
in  principle,  be  significantly  ameliorated  by  the  construction  of  lightweight  sandbar¬ 
like  structures,  which  could  be  "tuned”  to  the  most  damaging  waves,  thus  damping 
them  considerably.  This  technique  would  replace  the  present  heavy  and  very  ex¬ 
pensive  barrier  walls,  which  may  impinge  harmfully  on  the  natural  balance  in  the 
environment. 

•  The  above-mentioned  resonant  effect  may  also  be  used  to  produce  the  opposite 
effect:  the  bar-like  structures  could  act  as  a  lens,  concentrating  the  power  of  the 
most  energetic  waves  into  a  small  region  and  thus  increasing  the  efficiency  of  water- 
wave  driven  electric  generators. 

•  The  approach  taken  to  understand  the  formation  and  evolution  of  sand  ridges  may 
be  applicable  in  some  degree  to  other  structures  in  nature  that  are  the  result  of 
fundamentally  nonlinear  interactions,  such  as  cloud  patterns  and  sand  dunes. 

1.3  Morphology  of  Oceanic  Sedimentary  Structures 

Not  that  long  ago  it  was  thought  that  sand  ripples,  like  those  found  in  the  beach  zone,  and 
their  larger  cousins  the  sand  ridges,  were  morphologically  similar.  We  now  recognize  a 


variety  of  different  sedimentary  structures,  defining  the  categories  by  shape  or  generating 
mechanism.  Examples  are  sand  ripples,  ridge-runnel  systems,  tidal  ridges,  longshore 
sand  ridges.  The  formation  and  maintenance  of  these  sedimentary  structures  is  not  well 
understood. 

In  the  near-beach  zone,  including  the  breaker  zone,  occur  small  sand  ripples,  on 
the  order  of  a  few  centimeters  high,  which  come  in  a  multitude  of  shapes  and  forms. 
Larger  structures,  such  as  crescentic  bars,  occur  as  well.  In  this  region  the  fluid  flow  is 
quite  complex,  since  there  are  both  incident  and  reflected  waves,  tidal  flow  effects,  and 
turbulence  from  wave  breaking. 

The  ridge-runnel  system,  so  common  in  the  near-beach  zones  in  the  American  North¬ 
east  and  in  the  Great  Lakes  [4],  is  comprised  of  a  large  bump  3  to  15  meters  away  from 
the  beach,  about  0.3  meters  high  and  up  to  perhaps  7  to  10  meters  in  length,  which  is 
preceeded  by  a  runnel.  The  runnel  may  or  may  not  be  scoured  with  small  ripples.  The 
system  is  thought  to  be  formed  by  storms  eroding  the  beach  and  the  dune  fields  and/or 
by  tidal  currents  [4].  Davis  et  al.  [4]  provide  observational  evidence  for  their  claim  that 
storms  seem  to  play  a  minor  role  in  the  evolution  of  these  structures  once  they  have 
formed. 

Tidal  ridges,  which  were  noticed  by  Off  [5],  are  rhythmic  features  oriented  parallel 
to  the  direction  of  tidal  currents.  They  are  8  to  30  meters  high,  7  to  60  kilometers 
long,  and  separated  by  1  to  10  kilometers.  Allen  [6]  found  that  their  height  is  roughly 
proportional  to  the  square  root  of  their  spacing,  and  that  they  are  composed  of  sand,  silt, 
and  mud.  He  reported  that  they  occur  where  tidal  currents  reach  at  least  1  to  5  knots 
and  where  there  is  an  ample  supply  of  sediment.  Tidal  ridges  are  also  very  prominent  in 
the  neighborhood  of  river  deltas.  Tidal  ridges  may  have  a  fairly  flat  dome,  suggesting  to 
some  researchers  that  erosion  effects  play  a  very  minor  role. 
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Sandbars  are  distributed  in  complicated  patterns  on  the  continental  shelf,  and  it  is 
sometimes  difficult  to  discern  which  is  a  tidal  ridge  and  which  is  a  longshore  sand  ridge, 
the  object  of  attention  in  this  study.  For  example.  Figure  1.1,  taken  from  a  paper  by 
Swift  [7],  shows  the  relative  orientation  of  different  types  of  ridges.  Note  that  some  bars 
fan  out  around  river  deltas,  while  some  are  oriented  parallel  or  almost  normal  to  the 
coast. 


Figure  1.1:  Submerged  ridge  field  from  Long  Island  to  Florida,  from  Swift  [7]. 

Longshore  sand  ridges  are  common  features  of  the  continental  shelf  in  water  deeper 
than  the  surf  zone,  from  the  near-shore  region  to  the  farthest  reaches  of  the  shelf.  The 
better-known  ridge  fields  are  those  found  in  the  shallowest  end  of  this  range,  primarily 
because  they  are  readily  seen,  as  is  illustrated  in  Figure  1.2,  which  shows  the  bar  system 
off  central  Harrison  County,  Mississippi.  Other  near-shore  systems  are  found  along  the 


coasts  of  the  Carolinas.  Florida,  the  northern  coast  of  Alaska,  in  the  Black  Sea.  the  Baltic 


Sea.  and  even  in  large  lakes  such  as  Lake  Michigan.  Longshore  sand  ridges  can  be  found 
in  the  farthest  reaches  of  the  shelf  hugging  every  continent  around  the  world  as  well. 

From  observations  there  seems  to  be  a  mean  gradient,  in  the  neighborhood  of  0.02  to 
0.05.  which  favors  the  formation  of  longshore  sand  ridges  [8].  The  ridges  are  composed 
mostly  of  fine  sand  and  silt,  sometimes  of  mud.  The  mean  sediment  particle  size  ranges 
between  0.1  and  0.5  millimeters.  As  shown  in  Figure  1.3,  the  ridges  are  typically  a  few 


Figure  1.2:  Sand  ridges  in  shallow  water,  Harrison  County,  Mississippi. 

meters  high  and  are  spaced  hundreds  of  meters  from  each  other.  Groups  of  up  to  12 
ridges  have  been  seen,  mostly  parallel  to  each  other.  Their  migration  rates  vary  from 
place  to  place;  for  instance,  the  ridges  on  Sable  Island  Bank  have  been  estimated  to 
move  at  rates  ranging  from  0.5  meters  per  year,  in  water  60  meters  deep,  to  5  meters 


per  year,  in  30-meter  depths  [9].  The  ridge  fields  are  routinely  found  in  regions  where 
the  water  depth  is  small  compared  to  the  wavelength  of  surface  waves  with  frequencies 
in  the  infra-gravity  range  [10]. 
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Figure  1.3:  Cross  Section.  Sand  ridges  off  the  coast  of  Northern  Alaska.  Almost  half  of 
the  1350  Km.  long  coast  share  such  morphology.  From  Short  [10]. 

1.4  Comments  on  Field  and  Laboratory  Observations 

In  addition  to  the  inherent  difficulties  of  conducting  laboratory  experiments  involving 
liquid/sediment  media  (such  as  leveling  the  sediment  bed  after  each  trial,  extracting  gas 
bubbles  and  contaminants,  etc.),  laboratory  experiments  that  purport  to  mode!  oceanic 
phenomena  are  difficult  to  interpret  since,  in  most  cases,  oceanic  phenomena  do  not 
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scale  well  to  laboratory  conditions.  Field  observations  are  also  quite  challenging.  The 
environment  can  be  quite  hostile,  time  scales  are  long,  and  spatial  scales  are  large. 

It  was  not  until  the  mid- 1940s  that  exploration  into  the  deeper  parts  of  the  conti¬ 
nental  shelf  was  even  considered.  In  the  '60s  and  '70s  a  great  deal  of  field  observations 
were  made  on  sand  structures  of  ail  sorts.  Nevertheless,  ridge  fields  have  just  begun  to 
be  investigated  in  a  systematic  way.  Our  expertise  with  signal  processing,  telemetry,  and 
acoustical  and  radio  wave  remote  sensing  have  only  recently  been  upgraded  to  the  point 
where  large  scale  or  long  term  experiments  are  now  possible  [11].  Although  acoustical 
wave  remote  sensing  has  been  shown  to  be  the  best  way  to  probe  the  ocean  environment, 
we  lack  the  concerted  effort  that  would  be  required  to  produce  large-scale  acoustic  array 
measurements  that  would  enable  time-dependent  data  gathering  of  the  bottom  topogra¬ 
phy  and  fluid  motion.  It  is  no  surprise,  then,  that  very  few  complete  data  sets  of  ridge 
fields  exist  in  the  open  literature  at  present.  In  addition,  what  is  meant  by  a  •‘complete" 
data  set  has  cen  changing  over  the  years.  In  our  study,  a  complete  data  set  would 
include  bathimetric  records,  as  well  as  surface  or  internal  wave  directional  spectra  taken 
over  the  course  of  years. 

What  do  experiments  suggest  about  sedimentary  transport  in  an  oceanic  environ¬ 
ment?  In  the  case  of  laboratory  experiments  with  sand  ripples,  theory  seems  to  qualita¬ 
tively  agree  with  experiments  for  a  rather  limited  regime  of  flow  and  time  spans.  Some  of 
the  most  carefully  conducted  sand-ripple  experiments  are  those  of  Boczar-Karakiewicz, 
Benjamin,  and  Pritchard  [12].  However,  Pritchard  [13]  has  stated  that,  based  on  his 
as-ye*  unpublished  results,  with  an  erodible  bed  in  a  standing  wave  tank,  a  wave  field 
can  show  very  long  periods  of  homogeneous  activity  with  little  discernible  movement  of 
the  bottom.  Then,  at  an  unpredictable  moment,  if  all  is  right,  the  waves  can  grow  to 
the  point  of  breaking.  A  great  deal  of  turbulence  is  seen  in  the  boundary  layer,  the  bed 
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suffers  a  very  quick  rearrangement,  and  the  wave  field  returns  to  little  activity.  Pritchard 
did  not  measure  all  of  the  fluid  parameters  in  the  water  column  or  in  the  suspended  sed¬ 
iment.  In  his  opinion,  insofar  as  sand  ripples  in  the  near-shore  zone  are  concerned,  wave 
breaking  is  an  extremely  important  source  of  sediment  structure  formation. 

Is  wave  breaking  essential  to  the  generation  of  sand  ridges?  Ripples  can  be  formed  in 
a  laboratory  tank  by  a  non-breaking  wave  field.  Sand  ridges,  as  was  mentioned  earlier, 
appear  in  regions  where  no  breaking  waves  occur.  Most  models  for  the  near-  or  far- 
shore  zones,  like  the  one  which  will  be  presented  here,  do  not  apply  to  the  breaker  zone. 
While  breaking  is  an  excellent  source  of  turbulence,  we  do  not  know  how  it  controls  the 
dynamics  of  sediment  and  of  the  underlying  sand  structures.  Nevertheless,  wave  breaking 
away  from  the  breaker  region  has  been  seen  to  have  the  following  effects:  Lau  and  Travis 
[8]  found  that  sandbars  beyond  the  breaker  zone  do  not  disappear,  but  simply  change 
location  after  a  severe  storm.  Short,  in  his  field  observations  in  Northern  Alaska  [10], 
found  that  severe  storms  seem  to  rework  the  bars,  but  that  some  sandbars  photographed 
in  1949  and  1955  were  still  identifiable  after  approximately  30  years.  Preliminary  data 
from  the  so-called  “Super  Duck”  [11]  experiments  (purported  to  be  the  most  conclusive 
measurement  enterprise)  show  this  bar  “reworking”;  we  are  waiting  for  the  release  of 
these  data. 

There  are  two  main  differences  in  the  near  and  far  ends  of  the  continental  shelf  insofar 
as  the  fluid  environment  is  concerned.  First,  in  the  near-shore  we  can  identify  strong 
incident  and  reflected  components  to  the  wave  field.  Second,  as  the  (nonlinear)  waves 
shoal  some  of  the  energy  in  the  lower  frequencies  will  shift  to  higher  frequencies.  Not 
only  is  there  significant  asymmetry  in  the  velocity  field,  there  can  be  quite  pronounced 
asymmetry  in  the  acceleration  field.  Bijker  et  al.  [14]  made  laboratory  measurements  of 
acceleration  and  velocity  fields  for  water  waves  with  fairly  high  Stokes  numbers,  in  the 
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order  of  12-57.  They  found  the  net  transport  to  be  in  the  direction  of  the  wave,  particu¬ 
larly  if  the  wave  was  very  nonlinear.  Smaller  particles  seemed  to  be  transported  mostly 
by  the  Stokes  flow,  whereas  larger  particles  responded  mostly  to  the  "acceleration"  field. 
Hallermeier  [15]  analysed  a  large  experimental  data  set  and  found  an  empirical  rule  for 
the  prediction  of  ripple  characteristics  based  on  the  acceleration  field,  which  suggests  that 
this  field  may  be  an  important  sand-transport  mechanism  in  the  near-beach  zone.  Elgar 
et  al.  [16]  made  measurements  in  the  shoaling  region,  in  water  depths  in  the  range  of  1-6 
meters,  over  a  topography  with  mean  slope  of  5  %.  which  confirmed  the  existence  of  the 
velocity  and  acceleration  field  asymmetry.  They  found  that  the  acceleration  asymmetry 
becomes  increasingly  significant  with  decreasing  water  depth.  The  above  investigations 
suggest  that  the  acceleration  field  becomes  ever  more  important  as  the  distance  to  the 
beach  decreases;  our  model  would  not  apply  in  this  area,  since  the  transport  equation 
we  use  does  not  include  acceleration  effects. 


1.5  Sedimentation  Transport  Models 

As  mentioned  previously,  much  of  the  work  on  sedimentation  has  been  designed  to  under¬ 
stand  how  the  sediment  moves,  rather  than  how  it  generates  patterns.  Most  researchers 
working  on  sedimentation  transport  assume  an  outer  fluid  flow  at  the  edge  of  a  bound¬ 
ary  layer,  and  attempt  to  model  sediment  motion  on  the  bed  and  in  the  layer.  Sleath 
presents  a  good  review  of  the  subject;  we  will  summarize,  therefore,  only  in  a  cursory 
manner  the  different  sedimentation  models. 

A  model  developed  by  Bagnold  [17,18]  assumed  that  wave-induced  oscillatory  water 
motion  causes  sediment  to  move  back  and  forth  with  a  net  expenditure  of  energy.  Al¬ 
though  no  net  transport  results  in  such  an  oscillatory  flow,  the  energy  dissipation  acts  to 
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keep  the  sediment  in  suspension.  Once  in  suspension,  any  steady  current  superimposed 
on  this  oscillatory  flow  will  then  cause  a  net  transport  of  the  suspended  sediment  in  the 
direction  of  the  instantaneous  total  bottom  stress.  Originally  a  bed  load  model.  Bag- 
noid's  model  is  also  applied  to  suspended  load  transport  for  low  Froude  number  flows. 
A  threshold  of  motion  parameter,  called  the  Shield's  parameter,  is  incorporated  into  the 
model  to  reflect  the  fact  that  a  critical  amount  of  energy  must  be  imparted  on  the  bed 
before  transport  can  occur.  Smith  [19]  and  Fredsoe  [20]  applied  this  model  to  the  ocean 
environment.  They  assumed  a  constant  eddy  viscosity  and  obtained  criteria  for  the  on¬ 
set  of  instability  and  ripple  formation.  Richards  [21]  used  instead  a  turbulent  scale  that 
increases  linearly  in  height  from  the  bed.  thus  obtaining  two  modes  of  instability,  which 
yield  small-  and  large-scale  ripples  respectively.  Bagnold's  model  has  also  been  used  with 
some  success  in  the  near-shore  zone,  in  a  version  which  includes  the  effect  of  wind  on 
sediment  transport  rate  [22].  However.  Bailard  and  Inman  [18]  found  that  the  model  did 
not  perform  adequately  when  the  waves  are  not  normally  incident  to  the  beach. 

Another  sedimentation  model  by  Raudkivi  [23].  and  by  Williams  and  Kemp  [24]. 
attributes  the  formation  of  ripples  to  a  chance  piling  of  sediment.  This  deformation 
then  causes  the  flow  to  separate,  with  subsequent  building  up  of  the  ripple  downstream. 
They  attribute  the  initial  small  deformation  to  the  random  action  of  highly  turbulent 
velocities,  or  “bursts”,  close  to  the  bed. 

Lastly  we  mention  the  model  in  Longuett- Higgins'  seminal  paper  [25].  He  shows 
how  a  second  order  drift  velocity,  which  was  first  noted  by  Stokes  [26].  develops  in  the 
boundary  layer  from  an  outer  oscillatory  flow  or  in  the  bulk  of  the  fluid  through  the  action 
of  nonlinear  waves.  This  drift  velocity  is  capable  of  transporting  sediment,  particularly 
suspended  sediment.  A  number  of  people  have  studied  this  mechanism;  of  note  are 
Johns  [27],  who  developed  explicit  expressions  appropriate  for  the  ocean  environment 
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and  studied  the  character  of  .he  drift  velocity  and  its  stability,  and  Blondeaux  [28]  and 
Yittori  and  Blondeaux  [29],  who  looked  at  the  stability  and  formation  for  Froude  numbers 
at  which  flow  separation  does  not  occur.  They  determined  adequate  height,  spacing,  and 
onset  thresholds,  as  compared  to  laboratory  experiments.  The  second  of  these  papers 
introduced  more  structure  and  made  a  case  for  the  inclusion  of  nonlinear  effects. 

In  our  study  we  adopt  this  last  model.  The  mean  slopes  in  those  regions  of  principal 
interest  here  are  very  low,  hence  down-slope  gravitational  transport,  which  is  important 
in  the  coastal  environment,  plays  a  negligible  role  in  this  model.  The  ratio  of  bar  height 
to  separation  is  very  much  below  the  critical  value  of  0.1.  As  noted  by  Sleath  [3],  values 
above  0.1  usually  lead  to  boundary  layer  separation  behind  the  crests  of  the  bars,  and 
vortex  formation  takes  place.  When  this  occurs  vortex  ripples  will  spread  over  the  entire 
bed. 

1.6  Sedimentary  Bar  Models 

Among  the  researchers  who  have  coupled  a  sedimentation  transport  model  to  an  oceanic 
w'ave  field  to  look  at  the  process  of  bar  formation  in  the  oceanic  environment  are  Holman 
and  Bowen  [30].  They  use  the  linear  three-dimensional  water  wave  equations  to  compute 
drift  velocity,  which  in  turn  they  substitute  in  Bagnold’s  transport  model  for  suspended 
load.  In  particular,  they  examine  the  edge  wave  case  in  an  effort  to  compute  the  formation 
of  crescentic  bars  in  the  shoaling  region.  Bowen  [31]  has  also  examined  the  performance 
of  his  model  in  predicting  the  spacing  of  longshore  ridges  and  reports  good  qualitative 
agreement  with  field  observations. 

As  mentioned  earlier  in  connection  with  Pritchard's  work  [12],  laboratory  and  field 
observations  indicate  that  standing  wave  patterns  display  a  Bragg  resonance  process 
with  an  underlying  wavy  bottom.  In  the  steady-state,  the  ripples  develop  a  spacing 
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that  is  roughly  half  the  local  average  length  of  the  water  waves.  This  theory  [12.32-34] 
is  applicable  in  the  near-shore  environment,  since  it  relies  on  the  scouring  effect  of  a 
standing  wave  pattern.  This  first  order  theory  is  the  one  most  widely  studied,  since  it  is 
most  easily  implemented  in  the  laboratory:  at  one  or  another  time,  researchers  implicated 
this  mechanism  in  the  generation  of  all  sandbars. 

The  ridge  and  runnel  system  has  been  modeled  using  a  variant  of  Bagnold's  transport 
formula  by  Dean  [22]  and  deYriend  [35].  The  extent  of  their  success,  however,  is  hard  to 
discern  from  their  publications.  Since  the  undertow  and  the  local  bed  slope  are  significant 
and  since  the  effect  of  the  wind  in  generating  stresses  on  the  surface  of  the  ocean  must 
be  taken  into  account,  modeling  the  formation  of  runnels  is  difficult.  Russell  and  Osorio 
[36],  Bijker  et  al.  [37]  found  that  on  a  sloping  beach,  the  mass  transport  velocity  near 
the  bed  was  onshore  before  breaking  and  offshore  after.  This  effect,  independent  of  wave 
reflection  from  the  beach,  may  explain  why  these  bar  systems  are  usually  found  close  to 
the  plunge  line  of  breakers. 

Huthnance  [38]  develops  a  theory  for  the  formation  of  tidal  ridges,  based  on  an 
instability  which  is  triggered  by  a  small  protuberance  on  the  shelf.  The  ensuing  boundary 
layer  develops  a  bar  that  is  fed  by  bedload.  The  resulting  steady-state  bar  is  finite  in 
extent  and  parallel  to  the  always  present  currents.  Equilibrium  is  reached  when  the 
supply  of  sand  is  exhausted.  Huthnance  notes  that  the  tops  of  these  ridges  are  fiat  rather 
than  rounded,  which  he  claims  dismisses  erosion  as  being  the  source  for  the  generation 
of  these  structures.  Erosion  should  not  be  dismissed,  however,  since  these  bars  appear 
close  to  river  deltas  and,  possibly,  as  features  of  older  beaches.  Huthnance's  study  does 
not  address  the  periodic  nature  of  these  bars  and  does  not  suggest  a  relation  between 
their  height  and  spacing. 

Theories  for  the  formation  of  the  longshore  sand  ridges,  which  are  the  subject  of  this 
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study,  will  be  briefly  reviewed  in  the  following  section. 

1.7  Proposed  Model 

1.7.1  Historical  Development 

Among  the  first  to  suggest  that  infragravity  standing  waves  may  be  responsible  for  sand 
ridge  formation  was  Suhaida  [39j.  He  did  so  at  a  time  when  few  people  saw  anything 
fundamentally  different  about  near-shore  sandbars,  where  a  strong  standing  wave  field  is 
present,  and  bars  or  ridges  far  from  the  beach,  where  little  or  no  standing  wave  pattern 
is  to  be  found.  Short  [10]  made  field  measurements  of  sand  ridges  in  Alaska.  He  found  a 
loose  correlation  between  the  ridge  spacing  and  the  average  peak  infragravity  component 
wavelenth. 

Lau  and  Travis  [8]  derived  a  drift  velocity  from  a  Stokes  water  wave  field  for  a  bed 
with  constant  slope.  They  were  able  to  estimate  the  spacing  and  the  number  of  ridges 
from  the  periodicity  of  the  drift  velocity.  They  made  use  of  the  SRIT  (slightly  reso¬ 
nant  interacting  triads)  approximation  developed  by  Lau  and  Barcilon  [40]  and  Mei  and 
Uinluata  [32]  for  weakly  nonlinear  shallow  water  waves  to  solve  approximately  for  the 
wave  motion.  They  made  some  comparisons  with  field  data  to  examine  the  adequacy  of 
their  theory  in  predicting  the  observed  bar  separation  distances. 

Boczar-Karakiewicz  brought  this  problem  to  the  attention  of  Bona  while  the  latter 
was  visiting  Poland  in  the  early  1980s.  Eventually,  their  collaboration  resulted  in  the 
Boczar-Karakiewicz,  Bona,  Cohen  paper  [1],  in  which  they  use  the  ideas  of  Lau  and  Bar¬ 
cilon  to  obtain  the  resulting  drift  velocity  in  a  boundary  layer  and  use  this  drift  velocity 
as  a  source  of  sediment  motion  in  a  transport  equation.  Exploiting  the  discrepancy  of 
the  time  scales  between  fluid  and  sediment  dynamics,  they  were  able  to  formulate  the 
first  truly  evolutionary  sand  ridge  models.  Their  model  is  appropriate  for  the  shallower 
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end  of  the  continental  shelf,  since  it  was  derived  for  an  isotropic  water  environment. 
Later,  the  model  was  extended  to  the  internal  wave  case  and  was  tested  against  actual 
field  data  [2].  Encouraged  by  the  results  of  the  field  data  comparisons,  it  was  thought 
that  the  natural  extension  of  this  ongoing  research  project  should  be  to  increase  the 
model's  applicability  to  three  dimensions.  The  result  is  the  present  study.  Bona  and 
this  author  are  currently  pursuing  some  of  the  more  theoretical  iosues  in  the  project, 
while  Boczar-Karakiewicz  is  testing  the  models  against  field  data  and  is  investigating 
possible  practical  improvements  to  the  model,  such  as  the  use  of  more  realistic  transport 
equations  and  the  addition  of  more  phenomenology,  so  that  the  model  might  prove  useful 
to  the  engineering  community. 

1.7.2  Brief  Description  of  the  Model 

Referring  to  Figure  1.4,  we  envision  infra-gravity  waves  coming  into  the  purview  of  the 
model  at  the  line  x  =  0,  which  is  set.  on  the  deep  side,  by  the  point  at  which  the  long 
waves  “feel”  the  bottom  topography.  The  shoreward  direction,  x ,  increases  as  the  wave 
travels  shoreward.  The  span-wise  direction,  given  by  y,  is  approximately  parallel  to  the 
line  of  constant  phase  of  the  incoming  waves.  The  waves  propagate  shoreward,  possibly 
at  an  angle  with  respect  to  the  prevailing  direction  of  maximum  gradient  of  the  bottom 
topography.  In  the  deeper  reaches  of  the  shelf,  the  waves  would  be  supported  by  the 
picnocline,  while  in  the  isotropic  water  column,  the  waves  would  be  on  the  ocean  surface. 
The  extent  of  the  model  is  limited  in  the  longshore  direction  by  the  disintegration  of  the 
interface  supporting  the  internal  waves,  by  the  approach  to  the  breaking  zone,  bv  any 
singularity  in  the  depth,  and  by  significant  energy  transfer  from  low  to  high  frequencies. 
The  span-wise  direction  is  limited  by  the  same  sort  of  issues.  Taking  advantage  of 
the  disparate  time  scales  for  bottom  and  fluid  evolution,  the  assumed  gently  sloping 
bottom  will  be  considered  to  appear  as  a  fixed  but  non-uniform  surface  to  the  waves 
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as  they  progress  and  eventually  dissipate  on  the  shore.  This  assumption  enables  us  to 
decouple  the  problem:  starting  with  some  initial  bottom  configuration,  we  solve  the 
hydrodynamics  that  evolve  in  time  t  and  find  the  drift  velocity  in  the  boundary  layer; 
the  resultant  drift  velocity  is  then  used  in  a  transport  equation  to  update  the  bottom 
topography,  which  is  evolving  in  time  scale  T.  which  is  considerably  longer  than  t. 


Figure  1.4:  Aereal  view  of  the  problem. 

1.7.3  General  Comments 

A  few  comments  must  be  made  as  far  as  the  general  mechanism  for  longshore  ripple 
and  sand  ridge  formation  is  concerned.  If  a  standing  wave  pattern  exists  in  the  surface 
waves,  Linear  or  nonlinear,  the  scouring  effect  of  the  waves  generates  ripples  obeying  a 
Bragg  scattering  mechanism.  This  is  a  first-order  phenomenon.  Its  ability  to  influence 
the  shape  of  the  bottom  topography  relies  on  the  existence  of  both  a  reflected  and  an 
incident  wave.  As  we  move  further  seaward,  the  reflected  component  may  become  weaker 
and  weaker.  Yet.  we  find  large-scale  bars.  In  this  region  it  is  suggested  that  the  Bragg 
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mechanism  gives  way  to  the  second-order,  strictly  nonlinear  theory  that  we  present  in 
this  study.  Thus,  we  envision  that  both  mechanisms  operate  along  the  continental  shelf, 
but  in  the  near-shore  the  first  order  theory  is  prevalent,  while  in  the  deeper  reaches  the 
second-order  theory  prevails. 

The  second  order  theory  does  not  have  to  be  strictly  unidirectional.  However,  for 
very  mild  slopes  and  large  distances  from  the  shore,  the  reflected  component  is  pre¬ 
sumably  very  weak.  Hence,  if  the  bottom  is  not  restricted  in  this  way,  as  compared 
with  spatial  changes  in  the  surface  waves,  the  reflected  component  provides  a  great  deal 
more  structure.  For  the  surface  wave  case  the  reflected  component  is  still  significant 
close  to  the  shoaling  region  (which  is  the  extreme  end  of  the  purview  of  this  model); 
the  waves  are  assumed,  in  the  unidirectional  case,  to  dissipate  sufficiently  so  that  the 
reflected  component  is  negligible.  For  internal  waves,  the  issue  of  dissipation  is  relatively 
more  straightforward:  As  the  density  stratification  collapses  in  the  shallower  reaches  of 
the  shelf,  the  water  column  is  no  longer  able  to  support  an  internal  wave.  Incidentally. 
Boczar-Karakiewicz  et  al.  [2]  have  found  that  the  area  in  which  stratification  collapses 
agrees  with  the  limit  to  which  sand  ridge  fields  appear. 

The  frequency  range  in  this  model  is  limited  by  assumptions  of  shallow  water  wave 
theory,  i.e.,  long  wavelengths  compared  to  the  local  water  column  depth.  For  surface 
waves,  the  periods  range  from  0.5  minute  to  0.5  hour,  and  energies  in  the  order  of 
102  -  10s  J/m2.  For  internal  waves,  the  range  is  on  the  order  of  minutes  to  an  hour  in 
the  period,  and  energies  as  high  as  10 6 J/m2.  The  frequency  range  is  infra-gravitational 
for  both  the  surface  and  internal  waves. 

Internal  waves  in  the  above- mentioned  frequency  range  are  caused  by  such  things  as 
the  action  of  tidal  forces  on  the  stratified  fluid  flow  in  places  in  which  a  sudden  height 
change  in  the  bottom  topography  occurs,  such  as  that  at  the  edge  of  the  continental 
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shelf.  For  surface  waves,  on  the  other  hand,  the  generating  mechanism  is  less  obvious:  a 
distant  storm,  long-fetch  wind  effects,  or  tidal  forces.  No  provision  is  made  in  the  model 
for  external  forcing,  such  as  by  the  wind.  This  restriction  basically  limits  the  frequency 
range  of  the  surface  case  to  very  long  waves.  However,  some  of  the  longer  waves  are 
only  observed  in  the  farther  reaches  of  the  shelf,  in  areas  where  the  assumptions  of  an 
isotropic  water  column  are  hardly  realistic.  In  these  areas,  internal  waves  take  over. 

The  equation  which  models  the  surface  waves  in  this  study  is  a  highly  truncated  modal 
expansion  of  the  Boussinesq  system.  In  principle,  howc  r,  there  is  no  reason  why  the 
actual  Boussinesq  system  itself  could  not  be  used.  Elgar  and  his  collaborators  [41]  have 
examined  the  issue  of  the  recurrence  of  solutions  to  the  modally  truncated  Boussinesq 
Equation  numerically  in  the  Stokes  parameter  regime  of  0(1).  They  found  that  the 
two-mode  case,  which  will  be  used  in  this  study,  displays  recurrence-like  solutions  over 
a  great  many  wavelengths  of  distance.  As  the  number  of  modes  is  increased,  they  find 
that  the  recurrence  is  confined  to  fewer  and  fewer  cycles  the  more  modes  are  used.  In 
addition  they  find  that  initially  very  narrow  spectra  undergo  more  recurrence-like  cycles, 
before  the  spectra  flatten,  than  do  initially  broad-banded  spectra.  Their  conclusion  is 
that  recurrence-like  solutions  are  an  artifice  of  a  severely  truncated  modal  expansion  of 
the  Boussinesq  Equation. 

As  is  mentioned  in  their  paper,  other  researchers  have  studied  the  issue  of  recur¬ 
rence  of  solutions  in  such  equations  as  the  Nonlinear  Schrodinger  Equation  (NLSE) 
and  the  Zakharov  Equation  (ZE),  which  share  the  common  feature  with  the  Boussinesq 
Equation  that  they  all  undergo  0(  1)  energetic  transfers  between  their  modes  over  large 
times/distances.  Caponi  and  his  collaborators  [42]  found  numerically  that  in  the  ZE. 
depending  on  the  initial  conditions,  the  solutions  were  either  "periodic,  recurring,  tran¬ 
sitional,  or  chaotic."  In  connection  with  the  NLSE.  Weideman  et  al.  [43]  found  that 
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solutions  may  be  recurrent  or  chaotic,  depending  on  the  particular  structure  of  the  dis¬ 
cretization  used  in  its  numerical  solution  and  on  the  number  of  degrees  of  freedom.  If  the 
discretization  preserves  Hamiltonian  structure,  the  orbits  are  homoclinic.  Otherwise,  if 
Hamiltonian  structure  is  not  preserved,  for  a  few  degrees  of  freedom  the  discrete  NLSE 
behaves  entirely  differently  than  in  the  continuous  case.  As  the  number  of  degrees  of 
freedom  is  increased,  the  solutions  to  the  discrete  and  continuous  NLSE  converge,  as 
does  the  Hamiltonian  structure. 

In  conclusion,  then,  if  Elgar  and  his  co-workers'  findings  prove  correct  (we  are 
presently  addressing  this  issue  in  a  separate  study),  we  may  be  modeling  the  water 
waves  in  this  study  incorrectly.  However,  the  observations  of  Elgar  et  al.  do  not  weaken 
in  any  way  our  conjecture  that  weakly  nonlinear  shallow  water  waves  are  responsible  for 
the  formation  and  evolution  of  sand  ridges.  After  all,  there  is  more  than  ample  obser¬ 
vational  evidence  that  these  nonlinear  waves  travel  over  very  vast  spans  of  ocean,  i.e.. 
that  their  spectra  recurs  a  great  many  times  before  they  lose  their  coherent  shape,  over 
regions  where  sar  ’  ridges  are  a  prominent  feature  of  the  ocean  floor.  Certainly,  their 
findings  do  not  square  well  with  the  recurrence-like  solutions  that  internal  waves  are 
known  to  possess  over  very  large  spatial  scales.  Their  research,  if  verifiable,  leads  us  to 
conclude  either  that  there  is  something  inherently  wrong  with  the  modal  expansions  of 
the  Boussinesq  Equation  as  models  for  these  types  of  waves  or,  more  interestingly,  that 
their  findings,  along  with  those  of  Caponi’s  and  Weideman’s  groups,  are  pointing  out 
that  something  as  yet  not  understood  but  rather  fundamental,  is  awaiting  discovery  in 
connection  with  discretizations  of  nonlinear  evolution  equations  of  the  type  discussed. 
There  is  still  another  possibility:  It  could  be  that  the  loss  of  coherence  after  a  few  recur 
rent  cycles,  in  certain  situations,  is  responsible  for  the  interesting  fine  structure  observed 
in  actual  sand  ridge  fields. 


Chapter  2 

The  Hydrodynamics  of  the 
Water-Wave  Problem 

2.1  Preliminaries 

Owing  to  the  striking  similarity  of  the  typical  bar  spacing  to  the  length  scale  at  which 
energetic  interactions  among  the  most  significant  modes  of  shallow  water  waves  takes 
place,  we  believe  that  longshore  sand  ridges  are  formed  by  flows  in  the  boundary  layer 
which  are  generated  by  these  weakly- nonlinear  long  water  waves.  We  refer  to  these  waves 
as  “shallow  water  waves”  because  their  wavelength  is  considerably  greater  than  the  local 
depth  of  the  water  column  on  which  they  propagate.  An  appropriate  description  for 
these  waves  is  given  by  the  Boussinesq  System  [44]. 

In  this  section,  starting  from  mass  and  momentum  conservation,  the  equations  for 
long  weakly-nonlinear  water  waves  are  derived,  detailing  along  the  way  the  assumptions 
and  approximations  relevant  to  the  oceanic  environment. 

2.1.1  Conservation  of  Mass  and  Momentum 

Consider  a  function  p(r.t)  defined  in  a  time  dependent  set  Qt  C  ft3,  representing  the 
density  of  the  fluid  in  such  a  way  that  the  total  mass  of  the  fluid  body  m(Q o)  is  equal  to 
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and  is  constant.  If  this  invariance  holds  it  follows,  presuming  that  both  the  density  p 
and  the  velocity  u  are  €  that 

Dt  I  p(T.t)<Pr-  I  (Dtp  +  pV  ■  u)d3r,  (2.2) 

where  the  time  derivative  is  the  convective  derivative,  and  we  are  making  use  of  the 
divergence  theorem  and  Leibniz’  integration  rule.  If  this  invariance  applies  in  every 
subdomain  of  the  fluid  body,  then 

Dtp  +  pV  •  u  =  dtp+  V  •  (pu).  .  (2.3) 

When  it  is  assumed  that  the  density  of  the  fluid  element  does  not  change  (although 
it  may  different  for  different  fluid  elements),  the  above  simplifies  to 

V  •  u,  (2.4) 

and  we  say  that  the  fluid  is  “incompressible".  Note  that  incompressibility  is  not  a  prop¬ 
erty  of  the  fluid,  but  rather  of  the  motion.  It  amounts  to  assuming  that  the  volume  is 
preserved,  i.e.  its  flow  is  in  conditions  of  constant  volume.  Equation  (2.4)  applies  to  the 
case  considered  in  this  study. 

Momentum  is  conserved  as  well.  Conservation  of  linear  momentum  asserts  that 

Dl  k  pu,fiT  =  k  pS*T  +  Jaa,  * '  ad‘r '  (2'5) 

where  n  is  the  outward  normal  on  dQt.  f  encompasses  all  bodily  forces,  such  as  gravity, 
coriolis.  etc.,  and  the  contribution  of  contact  forces  enters  the  balance  through  the  stress 
tensor  N.  Expressed  in  another  way.  momentum  conservation  asserts  that 


u-  L  ^ = k  + k 


VKd3r. 


(2.6) 
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It  is  convenient  to  express  N tJ .  with  i  and  j  running  from  1  to  3  in  this  instance,  as  the 
sum  of  an  isotropic  part  - pbtJ ,  having  the  same  form  as  the  stress  tensor  in  a  fluid  at  rest, 
and  a  remaining  "deviatoric'’  part  dtJ  contributing  to  the  tangential  stresses.  The  tensor 
dtJ  has  the  distinctive  property  of  being  due  entirely  to  the  existence  of  the  motion  of  the 
fluid.  Furthermore,  the  deviatoric  tensor  may  be  recast  in  terms  of  physically  amenable 
terms.  Assume  that  dtJ  is  linearly  proportional  to  velocity  gradients,  so  that  the  stress 
tensor  is  now  =  -pbtJ  +  '2p(ttJ  3),  p  representing  the  pressure  (which  in  a  fluid 

in  motion  is  not  related  to  the  variables  of  state  in  equilibrium  thermodynamics),  and 
the  deviatoric  stress  tensor  separated  into  pure  straining  motion  and  expansion,  which 
are  respectively  the  second  and  third  terms  of  the  above  expression. 

In  almost  all  oceanic  circumstances  the  fluid  may  be  regarded  as  a  constant  density. 
Newtonian,  and  isotropic  fluid.  Thus.  =  eJt,  and  momentum  conservation  leads  to 
the  Navier-Stokes  equations 


where 


Du,  .dp  d  ,  .  ,  , 

p~dT  =  pf,~d7,+  d^{2p(etJ  ~ 


e  =i(^L  +  ^2) 

lJ  2  dxj  dx, 


(2.7) 


(2.8) 


is  the  rate-of-strain  tensor,  and  p  is  the  viscosity  of  the  fluid,  which  is  a  constant  of  pro¬ 
portionality  between  the  rate  of  shear  and  the  tangential  force  per  unit  area  when  >lane 
layers  of  fluid  slide  over  each  other.  The  viscosity  is  a  strictly  positive  quantity,  reflecting 
the  common  observation  that  the  force  between  layers  of  fluid  in  relative  sliding  motion  is 
always  a  frictional  force  resisting  the  relative  motion.  The  typical  value  of  p  for  water  at 
10°C  is  roughly  1.3T0 ~2gjcmj  ate.  and  it  decreases  at  about  3%  per  degree  centigrade 
rise  in  temperature  in  the  neighborhood  of  normal  temperatures.  When  appreciable 
temperature  differences  exist  in  the  flow  field,  p  must  be  considered  spatially  dependent. 
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However,  in  the  oceanic  setting  under  consideration,  such  temperature  differences  are 
not  present  and  the  viscosity  is  safely  assumed  constant.  Adopting  such  a  condition, 
and  incorporating  the  incompressible  condition,  conservation  of  linear  momentum  in  the 
bulk  of  the  ocean  is  expressed  as 


=  pf-  Vp  +  pV2u. 


(2.9) 


Viscous  terms  can  be  very  important  in  narrow  regions  of  flow,  or  in  very  small  scale 
motions,  where  the  significant  velocity  changes  are  confined  to  small  distances,  such  as  is 
the  case  in  the  so-calied  boundary  layers  at  the  air-water  interface  and  in  the  fluid  system 
immediately  above  the  bottom  topography.  Consideration  is  first  given,  however,  to  the 
effect  of  viscosity  away  from  both  boundary  layers  in  order  to  arrive  at  the  appropriate 
equations  of  conservation  in  the  bulk  of  the  oceanic  fluid. 

Given  that  it  is  the  fluid  is  isotropic  and  of  constant  density  p.  take  A,  the  wave¬ 
length  of  the  water  waves,  to  be  typical  of  the  length  of  appreciable  spatial  variation  in 
the  motion  or  magnitude  of  the  velocity  u.  Thus,  the  ratio  7 Z  =  pA|u|/p,  which  is  an 
appropriate  Reynolds  number  for  this  situation,  gives  an  estimate  of  the  relative  magni¬ 
tudes  of  the  inertial  forces  as  a  ratio  to  the  viscous  forces  involved.  For  long  wavelength 
waves,  { uj  may  be  replaced  by  the  more  accessible  velocity  measure  u;A.  where  &  is  the 
frequency  of  these  waves.  In  this  case  the  Reynolds  number  71  =  pa>A2/p  emphasizes  the 
fact  that  acceleration  in  the  fluid  is  proportional  to  frequency.  The  size  of  TZ  is  quite  iarge 
in  the  body  of  the  oceanic  fluid,  reflecting  the  fact  that  the  motion  is  almost  entirely 
governed  by  inertial  forces.  Therefore,  for  fluid  motion  which  is  dominated  by  inertial 
forces,  momentum  conservation  may  be  approximated  by 


=  pf-  Vp 


(2.10) 


to  reflect  such  balance  in  the  bulk  of  the  fluid.  Boundary  layers  form  on  the  air-water 
interface,  and  the  water-bottom  interface.  The  boundary  layer  at  the  water-bottom  in¬ 
terface  produces  significant  losses  even  in  the  ideal  situation  considered  here.  We  shall 
reserve  the  discussion  of  the  bottom  boundary  layer  for  a  later  chapter,  however.  The 
attenuation  due  to  dissipative  losses  in  the  air-water  interface  can  be  estimated  by  the 
following  argument:  Since  the  tangential  stresses  at  the  surface  are  zero,  and  the  normal 
stress  proportional  to  the  surface  tension,  the  losses  due  to  viscous  effects,  typified  by  the 
magnitude  of  v  =  p/p.  are  small  compared  to  the  inertial  forces,  that  is,  the  Reynolds 
number  1Z  =  *jk~2 /u  is  large  and  hence,  the  vortical  flow  may  be  neglected.  For  typical 
long  oceanic  waves,  1Z  ~  106  -  108.  Thus,  the  wave  decay  is  found  to  be  proportional  to 
ex p(~'2i>K2t)  from  the  conditions  imposed  on  the  tangential  and  normal  stresses  at  the 
surface,  where  n  is  the  wavenumber  of  the  waves.  For  the  typical  case,  the  "e- folding” 
distance  is  in  the  order  of  years.  This  simple  result  is,  of  course,  only  true  for  a  perfectly 
clean  and  wind-free  interface,  which  generally  is  not  the  situation  in  the  real  ocean  envi¬ 
ronment.  Wrhen  the  surface  is  contaminated,  the  free  surface  boundary  condition  is  more 
appropriately  modelled  by  an  elastic,  or  dynamic  no-slip  condition,  and  in  that  instance 
the  dissipation  is  not  trivial.  Additionally,  even  if  the  surface  was  clean,  the  energetic 
interactions  between  strong  wind  and  the  waves  usually  overwhelm  the  internal  frictional 
forces  just  discussed. 

The  contributions  of  coriolis.  surface  tension,  wind,  and  gravitational  forces  are  now 
briefly  examined.  With  the  exception  of  surface  tension  forces,  the  forces  just  mentioned 
enter  the  momentum  balance  as  terms  on  the  right-hand  side  of  Equation  (2.10).  The 
role  of  surface  tension  is  not  borne  out  of  conventional  momentum  balances  in  the  bulk 
of  the  fluid,  but  rather,  as  an  ad-hoc  condition  to  be  satisfied  at  the  air-water  interface. 
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The  apparent  body  forces  on  a  fluid  element  with  coordinates  at  rest  relative  to  a 
rotating  earth  with  approximately  constant  angular  velocity  il  are 

-2flxu-Qx(flxu).  (2.11) 

where  the  second  term  is  known  as  the  centripetal  force  and  the  first  term  as  the  coriolis 
force.  A  measure  of  the  relative  size  of  the  inertial  forces  to  the  coriolis  force  for  long  water 
waves  is  given  by  1/ LSI.  which  is  known  as  the  Rossbv  number.  Here  U  is  representative 
of  the  velocity,  L  is  a  characteristic  length,  and  Si  ss  7.3T0~;5sec~1 .  When  the  Rossbv 
number  is  greater  than  unity,  the  coriolis  force  is  negligible  as  compared  to  inertial  forces. 
For  long  waves  the  Rossby  number  may  be  estimated  by  the  ratio  uj/Sl.  and  this  ratio  is 
close  to  unity  for  planetary  waves.  Thus,  coriolis  forces  would  not  play  a  significant  role 
in  this  study,  since  the  wave  periods  for  the  waves  under  consideration  range  between 
fractions  of  a  minute  and  an  hour  -the  infra-gravitational  spectrum.  The  centripetal 
force  may  be  absorbed  into  the  pressure  term  in  the  Euler  equation,  i.e.  Equation  (2.10), 
or  safely  neglected  since  the  magnitude  of  the  inertial  forces  to  it  is  0(u>/Sl2). 

Continuity  of  stresses  at  the  air-water  interface  dictate  that  the  pressure  on  the  two 
sides  of  the  surface  can  differ  only  as  a  result  of  surface  tension  since  ideally  the  surface 
has  zero  mass.  The  force’s  origin  lies  in  the  fact  that  for  any  sufficiently  small  reversible 
isothermal  change  in  the  system,  the  total  work  done  is  proportional  to  the  gains  in  the 
total  Helmholtz  free  energy  [45].  The  molecular  origin  of  surface  tension  evidently  lies  in 
the  intermolecular  cohesive  forces.  The  magnitude  of  this  surface  tension  is  proportional 
to  the  local  interface  curvature.  Following  Batchelor  [45],  for  long  waves  the  pressure  at 
the  interface  is  regarded  as  constant.  The  condition  for  equilibrium,  with  z  the  height 
above  the  zero  level  reference  pressure  is 

1  1 

P9~  ~  ?(  +  s~)  =  constant, 

«i  ft  2 


(2.12) 
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where  Ri  and  are  the  radii  of  curvature  on  the  air  side.  The  constant  ~  is  the 
surface  tension  constant,  which  for  pure  water  at  15UC  is  about  lAdyn/ cm.  The  relevant 
parameter  which  reflects  the  relative  size  of  this  force  is  \J~]jpg.  For  pure  water,  the 
value  is  about  0.27cm,  which  is  the  length  scale  on  which  effects  of  surface  tension  are 
likely  to  be  comparable  to  the  effects  of  gravity.  Since  infra-gravitational  waves  have 
length  scales  many  times  larger  than  0.27cm,  surface  tension  effects  may  be  neglected. 

The  wind  is  a  major  source  of  energy  in  the  surface-wave  case.  The  water  wave 
spectrum  considered  in  this  study  spans  waves  with  periods  as  long  as  tens  of  minutes 
and  as  short  as  fractions  of  a  minute,  depending  on  the  water  depth.  The  longest  ones 
are  generated  by  earth  movements,  distant  storms  and  other  powerful  meteorological 
phenomena,  but  it  is  quite  clear  that  the  most  important  forcing  source,  for  waves  in  the 
high  end  of  the  frequency  range  of  the  model,  or  equivalently,  in  the  shallower  end  of  the 
shelf,  is  the  wind. 

Modelling  the  interaction  of  the  wind  with  the  ocean  surface  is.  at  present,  far  from 
satisfactory.  Nevertheless,  two  complementary  mechanisms  for  the  generation  and  main¬ 
tenance  of  waves  by  the  wind  have  been  proposed.  Phillip's  resonance  model  [46]  pro¬ 
poses  that  if  the  pressure  fluctuations  of  the  wind  are  in  phase  with  the  surface  waves, 
a  resonant  interaction  is  expected  with  ensuing  wave  growth.  It  is  said  to  govern  the 
initial  stages  of  wave  generation,  and  it  ignores  the  interaction  between  the  surface-wave 
and  the  actual  air  flow.  i.e.  it  considers  the  direct  action  of  turbulent  fluctuations  in 
aerodynamic  pressure.  A  complementary  theory  is  the  shear-flow  theory  worked  out  by 
Miles  [47],  which  governs  the  "instability"  phase.  In  this  theory  there  is  a  transfer  of 
momentum  from  the  wind  to  the  water  wraves  through  Reynold  stresses  at  the  boundary 
layer  with  ensuing  wave  growth. 

In  this  study  it  shall  be  assumed  that  the  wave  field  was  created  outside  of  the 
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model’s  region  of  applicability.  This  is  admittedly  a  gross  simplification  in  the  surface- 
wave  version  of  the  model,  but  again,  at  this  exploratory  stage  of  the  model,  consideration 
of  the  effect  of  the  wind  takes  us  far  from  more  immediate  issues. 

2.2  Surface  Wave  Problem 

By  “hydrodynamic  problem"  we  shall  mean  the  problem  in  the  domain  that  excludes 
the  boundary  layer  that  hugs  the  bottom.  As  we  shall  explain  later,  we  shall  exploit 
the  tr^  aendous  discrepancy  between  the  time  scales  of  the  fluid  motion,  represented 
by  t ,  and  the  time  scales  of  the  bottom  evolution  given  by  T.  The  latter  time  scale 
is  to  be  considered  a  parameter  in  what  follows.  As  illustrated  in  Figure  2.1.  one  may 
define  the  domain  for  the  hydrodynamic  problem  as  =  R2  x  [ -H(T )  +  £,//(<)]  % 
R2  x  [-H(T),  /?(<)],  since  6  <  |  H  |.  The  fluid  is  subjected  solely  to  gravitational  forcing. 
The  velocity  field  is  given  by  (u.  in),  where  the  first  entry  is  the  transverse  velocity  and 
w  is  the  vertical  velocity.  Position  is  represented  by  the  vector  (r,  z).  The  free  surface 
is  given  by  z  —  r)(r,t)  and  the  bottom  by  r  =  -H(r,T).  We  adopt  the  convention 
throughout  this  study,  that  the  operator  V3  =  V  +  kdz. 

It  is  postulated  that  the  fluid  is  initially  irrotational.  That  is. 

V3  x  (u.tt')  =  0.  (2.13) 

The  curl  of  Equation  (2.10)  yields  what  is  commonly  referred  to  as  the  Helmholtz  equa¬ 
tion: 

a 

— V3  x  (u,u?)  +  V3  x  (V3  x  (u,  u?)  x  (u,  u;))  =  0,  (2.14) 

at 

making  use  of  the  fact  that  the  force  field  is  conservative.  Appealing  to  Equation  (2.13). 
and  using  the  vector  identity  V3  x  V3  x  (u,  w)  =  — (u.  w)  x  V3  x  (u,  w ),  it  is  seen  that 


— (V3  x  (u.  tc))  =  V3  x  (u.  ic)-V3(u.  u’). 


(2.15) 
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z  =  -H(x,T) 


Figure  2.1:  Side  view,  surface  wave  problem. 

Since  V3  x  ( u,  w)  =  0  is  a  possible  solution  of  this  equation,  it  follows  that  if  the  flow  is 
initially  irrotational,  it  shall  remain  as  such  for  all  time. 

Since  the  fluid  motion  is  irrotational.  the  velocity  may  be  expressed  as  the  gradient 
of  a  scalar  potential  ®: 

( u,  w)  =  V3<p.  (2.16) 

From  conservation  of  mass.  Equation  (2.4),  the  equation  of  motion  within  the  fluid  is 

=  0.  in  Q.  (2.17) 

At  the  air-water  interface,  conservation  of  momentum  requires  the  pressure  to  be 
continuous.  The  assumed  constant  value  of  the  pressure  immediately  above  the  water  is 
set  to  zero.  Hence,  the  pressure  at  the  interface  when  the  surface  is  quiescent  shall  be 
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zero.  This  dynamical  boundary  condition  then  specifically  states  that 

4>t  =  |V3<A|2  -  grj,  at  ^  =  77.  (2.18) 

The  bottom,  which  is  assumed  impermeable,  has  a  normal  velocity  that  agrees  with 
that  of  the  fluid.  Thus, 

<pz  =  —  V//-Vd>,  at  2  =  —  K.  (2.19) 

Lastly,  the  kinematic  condition  on  the  air-water  interface,  that  fluid  particles  on  the 
surface  shall  remain  there  for  all  time,  may  be  expressed  as 

<t>,  =  T)t  +  V^-V7?,  at  z  =  77.  (2.20) 

Equations  (2.17)-(2.20)  with  the  additional  requirement  that  \  V<f>\  —  0  as  |r|  —  oc 
comprise  the  hydrodynamic  problem. 

Wave-like  solutions  of  the  above  boundary  value  problem  can  easily  be  derived  if  the 
waves  are  infinitessimal  in  amplitude.  Solutions  of  the  form  exp{i(K  r  -  cj(K)t)}  exist 
for  the  linearized  version  of  the  system,  provided  that  the  relation  between  the  frequency 
u  and  the  wavenumber  K  is  the  “dispersion  relation” 

u2  =  pKtanh(Klf ),  (2.21) 

where  k  =  |K|,  and  it  is  understood  that  the  dispersion  is  spacially  dependent  since 
H  =  H(r,T).  Both  k  and  u  must  be  real  if  we  are  strictly  interested  in  plane  wave 


solutions. 
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2.2.1  Hamiltonian  Formulation  of  the  Hydrodynamic  Problem 

To  obtain  a  Hamiltonian  formulation  of  the  hydrodynamic  problem  proposed  above,  it 
is  noted  that  the  motion  of  the  entire  fluid  body  will  be  determined  once  the  free  surface 
motion  is  known.  Specifically,  if  the  function  r)  which  describes  the  free  surface,  and  the 
velocity  potential  at  the  free  surface, 

$  =  <t>(r,z  =  r),t),  (2.22) 

are  known,  then,  for  each  t,  r?  determines  the  domain  and  $  determines  the  cor¬ 
responding  <t>,  which  is  the  unique  solution  to  the  hydrodynamic  problem  comprised 
of  Equation  (2.17)  through  Equation  (2.20)  and  Equation  (2.22),  with  the  additional 
assumption  that  |V3d>j  — <•  0,  as  |r|  —  oo.  In  what  follows,  we  rely  heavily  on  ideas 
developed  by  Zakharov  [48],  Miles  [49],  Bowman1  [50]  and  especially  Benjamin2  [51,52]. 

Consider  the  Hamiltonian  E  =  E(t),$).  The  choice  of  label  E  reflects  the  fact  that 
the  Hamiltonian  for  this  problem  is  conserved  and  is  numerically  equal  to  the  sum  of  the 
potential  energy  V,  and  the  the  kinetic  energy  K.  As  shown  by  Benjamin  and  Olver  [52], 
the  requirement  that  E  be  stationary  with  respect  to  independent  variations  of  6$  and 
b t]  yields  the  following  Hamiltonian  system: 

6E 

m  ~  J* 

6E 

*<  =  (2-23) 

where  the  derivatives  are  Gateaux  derivatives.3 

’S.  Bowman  has  unfortunately  left  science  and  is  now  an  actuarian  in  England.  He  leaves  behind 
some  good  work,  some  of  it  unpublished. 

2 Not  to  be  out  done  by  particle  theorists,  Benjamin  has  been  -successfully  so  far-  working  on  unifying 
fluid  dynamics  under  a  consistent  mathematical  framework.  A  possible  name  for  his  framework  could  be 
Super  Hamiltonian  Unified  Theory. 

3The  first  variation  of  the  functional  F  in  the  direction  r  is  defined  here  by  -jgF{ r  +  «r)e=o  ~  r £F, 
where  £  is  the  Eulerian  derivative. 
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To  briefly  demonstrate  the  equivalence  of  the  Hamiltonian  system  and  the  hydrody¬ 
namic  problem  consider  the  following:  The  energy  in  the  system  is 

v 

E  =  T  +  V  =  J  d2r  J  I|V3d>|2d.+  j  d2rl-g  t?  (2.24) 

R3  -«  R 2 

exactly.  Extremizing  E,  keeping  $  constant  and  considering  variations  6r), 

6E=  J  d2rfrq{^  \^3<t>\l=v  +  grj},  (2.25) 

R3 

which  leads  to  Equation  (2.18).  Next,  keeping  77  constant  and  considering  variations  6$, 
using  Equation  (2.17),  and  applying  Green’s  theorem, 

/oo  ( 00 

^  dEs<5$(Es)nSj-V^|Sj  +  J  ^  d£B6$(£B)n£s-Vd>|£B,  (2.26) 

where  n  is  the  outward  normal  to  the  boundaries.  Since  d£,  =  y/\~+J^rf)1d2r,  and 
d£ b  =  x/TT(WFd2r,  the  boundary  contributions  to  Equation  (2.26)  are 

J  d2rS${<frz -Vii-V<t>)\z=T„  (2.27) 

R3 

which  leads  to  Equation  (2.20),  and 

J  d2rt${<fiz  +  VH-V<j>}\z=_H,  (2.28) 

R3 

which  implies  Equation  (2.19). 

2.2.2  Development  of  the  Boussinesq  System 

The  Hamiltonian  System,  Equation  (2.23),  shall  be  specialized  for  the  case  of  weakly 
nonlinear  shallow  water  waves.  Define  the  parameters  a  <  1,  and  (3  <  1,  their  precise 
meaning,  in  terms  of  physically  relevant  parameters,  shall  be  made  clear  subsequently. 
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Assume  that  0(a)  ~  0(f32),  and  take  H  =  0(1),  V H  =  0(a).  tj  -  0(a),  $  =  0(a). 
Further,  consider  the  differentiations  d2,dt.V  =  0(/3). 

An  approximation  to  4>,  which  satisfies  the  boundary  value  problem  is 
<t>(r,  z,t)  =  $(r,<)  -  ±z2V2$(r,f)  -  zV-(HV$(r,  t)) 

0(a)  0(a/32)  0(a/32) 

which  can  be  easily  derived  using  Rayleigh’s  trick4  [53].  The  gradient  of  the  above 
expression 

U(r,  z,  t)  =  u(r,f)  -  {zV[V-(^u(r,f))]+  ^>2V(V-u(r,0)}  (2.29) 

gives  the  velocity  anywhere  in  the  inviscid  domain  of  the  fluid. 

The  potential  energy  is  exactly 

v  =  /  d2r^gri2.  (2.30) 

R2 

The  kinetic  energy  will  be  calculated  using  the  approximation  developed  above  for  the 
velocity  potential  developed  above,  Equation  (2.29): 

K  =  J  d2r{l-(H  +  ^(V*)2  +  |(Vtf -V*)2  -  ^(V2$)2},  (2.31) 

R2 

which  is  an  expression  of  0(a3/92),  and  0(a2/?4). 

Thus,  in  terms  of  the  velocity  at  the  surface  u  =  V$,  and  the  displacement,  the 
energy  is 

E  —  V  +  Aq  +  a  A  i  +  •  •  • ,  (2.32) 


’The  trick  was  first  used  in  connection  with  the  solution  of  the  electrostatic  field  in  an  axisymmetric 
strip,  unbounded  in  x,  say,  and  bounded  by  smooth  but  spatially  varying  edges  in  y.  The  harmonic 
functions  9  and  ii/  representing  the  potential  and  the  stream  function  are  expanded  as  <p  =  cos(ydx)f. 
and  tl>  =  svn(ydx)f,  where  f  is  determined  termwise  in  the  ■  mansion  in  terms  of  the  boundary  conditions. 
Thanks  are  due  to  Prof.  T.  B.  Benjamin  for  showing  me  this  trick,  and  for  introducing  me  to  all  matters 
Hamiltonian. 
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where 

A'o  =  J  d2r^Hu2  (2.33) 

R2 

^  +  (2.34) 

R2 

and  V'  is  as  before.  Substituting  E  in  Equation  (2.23),  to  lowest  order,  yields  the  wave 
equation 


th  +  V-(tfu)  =  0  (2.35) 

ut  +  gVr)  =  0.  (2.36) 

To  the  next  order, 

th  +  V-[(H  +  /?)u]  4-  V  [u V(ff2)VJ?  4-  \V(H3V- u)]  =  0  (2.37) 

ut  4-  (u-V)u  4-  5V77  =  0,  (2.38) 

a  version  of  a  Boussinesq  System  [44].  The  Boussinesq  System  (BSS)  is  a  shallow  water- 
long-wavelength-  weakly  nonlinear  approximation  to  the  Euler  Equation  which  admits  bi¬ 
directional  waves  as  solutions.  The  version  given  by  Equation  (2.37)  and  Equation  (2.38). 
however,  has  a  couple  of  troublesome  characteristics  from  the  standpoint  of  modelling 
a  physical  situation.  Namely,  the  system  is  linearly  unstable,  and  secondly,  it  is  rather 
poor  in  conveying  accurately  the  full  dispersion  relation. 

The  first  problem  may  be  seen  as  follows:  without  loss  of  generality,  consider  the  one¬ 
dimensional  version  of  BBS,  dropping  all  nonlinear  terms,  setting  g  =  1,  and  considering, 
for  simplicity,  the  case  of  a  uniformly  flat  bottom.  Additionally,  let  u(k,  t)  be  the  Fourier 
transform  in  space  of  u(x,t).  Cross  differentiating  and  combining  Equation  (2.37)  and 


Equation  (2.38),  the  resulting  equation  is 


uxx 


~  Him  —  Oi 


(2.39) 


or  equivalently, 

««(«)  =  (~«2  +  ^«4)u(k),  (2.40) 

for  which  a  solution  is 

u(n,t)  =  u(k,  0){Aexp(\J^K4  -  K2t)  +  5exp(-y^K4  -  n2t)}.  (2.41) 

It  is  then  immediately  obvious  that  the  solution  can  grow  «  exp(y^«;2t).  The  second 
problem  is  that  the  dispersion  relation  satisfied  by  Equation  (2.39), 

w2  -  k2(1- ^k2)  =  0  (2.42) 

is  an  adequate  approximation  to  the  Equation  (2.21)  strictly  for  very  low  frequencies  as 
shall  be  demonstrated  graphically  in  a  subsequent  section. 

2.2.3  Regularization  and  Scaling 

An  ad-hoc  procedure  which  “regularizes”  BBS  shall  enable  us  to  proceed  in  our  devel¬ 
opment.  As  an  alternative  model  to  the  Korteweg-deVries  equation  (KdV),  Peregrine 
[54]  developed  an  equation  which  has  eventually  been  referred  to  as  the  “Regularized 
Long  Wave  Equation”.  While  Peregrine  was  the  first  to  propose  it  as  an  alternative  to 
the  KdV,  it  was  Benjamin,  Bona,  Mahony  [55]  who  later,  but  independently,  proposed  a 
regularized  version  of  KdV,  with  the  express  intention  of  overcoming  some  of  its  short¬ 
comings  in  modelling  water  waves,  and  who  studied  the  resulting  model’s  properties. 
The  trick  is  to  use  the  lowest  order  continuity  and  momentum  balances  in  the  higher 
order  dispersive  terms  in  order  to  obtain  an  equation  which  is  more  amenable  to  physical 
and  computational  studies.  This  procedure  is  justified  on  the  grounds  that  the  error  in 
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making  the  substitution  into  the  dispersive  term,  which  is  a  higher  order  term,  shall  be 
no  worse  than  the  error  already  present  in  the  original  system.  We  shall  employ  a  similar 
technique  here,  exploiting  the  specific  form  of  the  bottom  topography,  being  careful  not 
to  destroy  the  bi-directional  nature  of  the  wave  solutions. 

Using  Equation  (2.36),  and  the  fact  that  VH  —  0(e).  approximate 

V-[uV(ff2)VII  +  \v(H3V-u)}  =  -\v-[V(H\))  +  0(a,e).  (2.43) 

t>  ») 

Thus,  the  regularized  system  (RB)  adopted  in  this  study,  as  an  approximate  model  for 
the  water  waves,  is 


rjt  +  V-[(H  +  tj) u]  -  ±V.[V(tf  2»fc)]  =  0  (2.44) 

ut  +  (u-V)u  +  gVi)  —  0.  (2.45) 

Several  comments  are  in  order.  First  of  all,  the  dispersion  relation  of  RB  is 

=  (2'46) 

For  plane  wave  propagation,  both  uj  and  «  must  be  real.  A  comparison  of  Equation 
(2.42)  and  Equation  (2.46)  shows  that  the  upper  limit  for  w  real  is  |k|  =  \/3  for  the 
BSB  system,  while  there’s  no  limit  on  RB.  Incidentally,  another  alternative  Boussinesq 
System,  the  version  used  by  Lau  and  Barcilon  [40]  results  in  a  dispersion  relation 

w2  -  k2  4-  =  0.  (2.47) 

This  dispersion  not  only  has  a  cutoff,  but  is  also  not  symmetric  about  k  -  0,  which  implies 
incorrectly  that  wave  propagation  in  the  forward  and  backward  directions  are  different. 
Figure(2.2)  illustrates  how  the  dispersion  relations,  Equation  (2.42)  and  Equation  (2.46), 
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compare  to  the  full  gravity  wave  dispersion  relation.  Equation  (2.21).  Both  H  and  g  were 
set  to  unity  and  j3  =  0.5  in  the  figure. 


K 

Figure  2.2:  Comparison  of  the  full  water  wave  dispersion  relation: - ,  BSB:  - 

. ,  RB: - . 

Regularization,  as  performed  here,  does  bring  a  couple  of  subtleties  that  must  be 
kept  in  mind:  (1)  Since  we  are  using  the  lowest  order  relations,  we  are  assuming  that 
the  solutions  are  wave-like.  However,  there  is  no  reason  to  expect  that  the  solution  will 
approximate  a  traveling  wave  solution  for  arbitrary  Cauchy  data.  (2)  We  have  preserved 
two-way  wave  propagation.  In  most  instances  when  regularization  is  performed,  the 
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Boussinesq  equation  that  results  is  applicable  to  strictly  one-way  wave  motion.  (3)  Since 
the  velocity  is  in  terms  of  the  surface  values,  rather  than  in  terms  of  averaged-depth 
velocity,  say,  the  irrotational  condition,  with  u  and  v  being  respectively  the  shoreward 
and  span-wise  velocity  components,  remains  in  the  simple  form 

Uy  —  t-V,  (2.48) 

which  is  quite  convenient  in  the  development  of  three-dimensional  problems. 

In  order  to  scale  the  hydrodynamical  model  developed  in  the  previous  section,  define 
a  parameter  related  to  the  size  of  the  amplitude  of  the  disturbance,  and  another  related 
to  the  size  of  spatial-or  temporal-changes.  The  first  one  is  a  =  a/h0 ,  which  gives  an 
estimate  of  the  degree  of  nonlinearity  in  the  problem.  The  second  is  02  =  (A //i0)2«  which 
conveys  the  degree  in  which  dispersive  effects  are  important.  The  Stokes  number,  which 
is  a  measure  of  the  balance  between  nonlinear  to  dispersive  effects  is  defined  as  the  ratio 
a/l32.  For  V  <  1,  nonlinear  effects  are  weak,  and  only  a  small  portion  of  energy  transfer 
occurs  on  moderate  space-time  scales,  so  that  0(1)  nonlinear  effects  are  possible  only 
after  very  large  scales.  For  U  ~  1,  inertial  effects  are  of  the  same  order  as  dispersive 
effects. 

Using  the  convention  in  what  follows  that  new  <—  scale  x  old ,  the  scaling  adopted  is 

->->!/«  r-f  >2'W> 

where  ho  is  a  characteristic  depth  of  the  water  column. 

In  addition,  we  seek  to  scale  span-wise  dependence  to  reflect  the  fact  that  we  are 
interested  in  a  case  of  propagation  that  is  primarily  shoreward  directed.  To  do  so.  it  is 
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assumed  that  there  is  a  const  <C  1  such  that 

0(|  i-K  |)  =  const  x  0(  |  y-K  |),  (2.50) 

for  which  a  consistent  uniform  expansion  of  the  RB  exits,  and  that  is  physically  relevant. 
It  may  be  shown  that  the  size  of  the  constant  is  of  the  order  of  the  reciprocal  of  3.  Since 
this  parameter  has  considerable  nuisance  value,  the  parameter  shall  be  set  const  —1/3 
for  the  rest  of  this  study.  This  implies  that  the  span-wise  variables  must  be  scaled 

y  <—  QJ/2y  y-u  «—  a~l^y-u,  (2.51) 

which  shall  alter  the  regularized  system,  but  shall  not  affect  the  irrotational  condition, 
Equation  (2.48). 

2.2.4  Description  of  the  Bottom  Topography 

Field  data  from  the  continental  shelf  suggests  that  there  are  two  time  scales,  a  fast  time 
scale  t  which  measures  the  evolution  of  the  fluid  quantities,  and  a  large  time  scale  T 
which  measures  the  evolution  of  the  bottom  topography.  In  addition  the  data  suggests 
that  the  typical  height  and  slopes  of  the  longshore  sand  ridges  is  such  that  e  =  0(Xnh)  = 
0(a).  Furthermore,  the  type  of  longshore  sand  ridge  under  consideration  is  such  that 
the  measure  of  longshore  spatial  variation  is  larger  than  the  spatial  variations  of  the  fluid 
quantities.  It  is  proposed  that  the  sand  ridge  shoreward  variation  be  X  =  ax.  Hence, 
two  scales  of  shoreward  variation  exist,  so  that 

dx^dx  +  adX-  (2.52) 

Thus  the  bottom  in  scaled  variables  is 


h(X,y,T)=  1  +ef(X,y,T), 


(2.53) 
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where  the  function  /  =  0(1). 

2.2.5  Slightly  Resonant  Interacting  Triads 
By  substituting  a  uniform  expansion  of  the  form 

7?  =  /o  +  Q1  /l  +  q2/2  +  ■  •  • 

u  =  go  4-  c^gi  +  a2g2  +  •  •  •  (2.54) 

into  Equations(2.44),  (2.45),  and  (2.48),  matching  order  by  order  we  are  able  to  solve  for 
the  surface  quantities  to  lowest  orders  in  a. 

For  the  momentum  equation,  Equation  (2.45),  the  first  three  orders  are 

a0  :  uot  +  »7ox  =  0 

W0(  +  lay  =0 

a1  :  Ult  +  u0u0x  +  rjlx  +  7/ox  =  0 

(2.55) 

Vlt  +  UQV0x  +  7/ly  =  0 
a2:  U2t  +  UlUOx  +  U0Uix  +  UqUqX  +  V0U0y  +  T)2X  +  T)ix  =0 

V2f  +  UlWOi  +  UoWlr  +  UoWoX  +  7/2y  =0 

Similarly,  for  the  irrotational  equation,  Equation  (2.48),  we  have 


a0: 

Woy  W02- 

=  0 

a1  : 

Uly 

-  Vix  -  Vox 

=  0 

(2.56) 

a2  : 

U2y 

-  v2x  -  ViX 

=  0. 

Finally,  Equation  (2.44)  yields 

<*°  ■  Vot  +  «0r  _  QfVOxxt  —  0 

a1:  tju  +  tii*  ~  ^-Vixxt  =  Fi(Tio,uo,vo,G-,x,X,y,t) 

a2:  r/2t  +  u2i  -  4"%rxt  =  F2(iio,  «o-  w0,  t?i,  «i,  »i,  X,  y,  t), 


(2.57) 


41 


where  G(.Y,yJ)=M 


r-  ,  3  VOyyt  .  23  (l)OxXt  H"  GtJ Oxr()  ,0  co, 

Fl  =  -t’Oy  +  - 2^  -  Uqx  -  UoTfox  -  Gu0x  -  TfoUOx  +  - - -  (2.58) 


F2  =  -vly  +  -  3:^;  -  UlX  -  UlT?0l  -  Gulx  -  7?it/0r  +  2-3  {.Tll*X^  +  Gr,lI*‘) 


—  UQlhx  ~  f/O^lx  GyV 0  —  Gxtlo  ~  ^O^Oy  “  Gt^Oy  —  VovOy  + 

-u0rj0X  -  Guqx  -  Vouox  +  +  4j--r0lt  + 


2 P  (Gyot) 


(2.59) 


Making  use  of  the  appropriate  irrotational  condition  when  a  simplification  is  possible, 
cross-differentiating  the  momentum  and  continuity  equations,  and  combining  the  results 
into  a  single  equation  yields 


a0  :  Crjo  =  0 

a1:  £77!  =  Giim,  «o,  v0,  G\  x,  X,  y,t) 

a2:  Crfr  =  vo,»7i,Tii,Vi,G;z,X,j/,t), 


(2.60) 


where 


C  —  dtt  ~  dxx  T'dxxtt* 


(2.61) 


Note  that  £  is  a  linear  operator  that  shows  up  at  every  order.  The  inhomogeneous  terms, 
G\  and  Qi  are,  respectively, 


Q\  —  (1  +  32dtt/3)r)oyy  +  G(1  +  2@2dtt/3)r)oxx  4-  2(1  +  32dtt/3)iloxX 
+  («o/ 2)r  -  («0 no)rU 


(2.62) 


42 


Qi  —  ( 1  +  32dtt/^)rhyy  +  0(  1  +  2/32(?tt/3)77irx  +  2(  1  +  3idtt/3)V\iX 
( 1  4-  t32dtt/3)voxx  +  0(1  +  2/325(t/3)r;oVy  +  2(1  4-  f32dtt/ 3)VoxX 
+Gx(l  4-  4{32dtt/ 3)  +  t)ox  +  GyTjQy  4-  2f3^GyyT}0tt/Z  4-  432GyTjotyy/3  (2.63) 

-(7?1“0  +  VoUi)xt  +  (uiUo)xi  +  («o)rX  -  («0%).Yi  +  G{uH 2)xx 
4"(uo/2)yy  "I"  (vo/2)xr  4"  (^?o/2)yj/  —  (rJOtvo)y 

Beyond  this  order  the  calculation  of  the  uniform  expansion  is  invalid  since  it  is  beyond 
the  order  to  which  RB  is  an  approximation  to  Euler’s  Equations.  The  lowest-order 
theory  is  suitable  since  we  are  only  interested  in  phenomenology,  rather  than  engineering 
accuracy.  Appendix  A  gives  the  reader  an  idea  of  the  sort  of  things  to  expect  at  0(a1). 

The  order  of  RB  and  the  two-scale  technique  restrict  the  region  of  validity  of  the 
present  model.  The  lowest  order  solutions,  which  are  linear,  are  valid  for  distances  that 
are  less  than  0(1 /ka),  the  scale  over  which  triads  of  Fourier  modes  exchange  significant 
energy.  Higher  order  terms  and  processes  neglected  in  the  expansion  restrict  the  range 
of  the  present  nonlinear  solutions  to  distances  less  than  0(1 /ha2).  Thus,  RB  is  not 
formally  valid  for  very  long  evolution  distances.  Boussinesq  equations  are  strictly  valid 
for  U  =  0(1),  but  they  are  quite  robust  [41].  In  this  study,  the  value  of  U  is  in  the  range 
of  10  to  30. 


Assume  the  shoreward  velocity  is 


u(x,X,y,t) 


Ho[a(fc,  X,  y)  +  aA(k,  X,  y)]e'lk*-^dk  +  c.c. 

+  JTooW*’  X,  y)  4-  aB(k,  X,  y)}e'Gkr-u;t)dk  +  c  c 


(2.64) 


where  k  =  i-K,  and  further,  assume  that  such  solution  may  be  approximated  by  the 
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discrete  Riemann  sum 

u(x,X,y,t )  =  y)  +  aAj(X,  -f  c.c. 

(2.65) 

+  L'£i  [bJ(X,y)  +  aBJ(X,y)}c^x-^  +  c.c.. 
where  c.c.  stands  for  complex  conjugate  of  the  expression  immediately  preceeding  its 
appearance.  The  a's  are  the  complex  incident  wave  amplitudes,  and  the  b's  are  the 
complex  reflected  wave  amplitudes.  The  reality  of  the  physical  variables  implies  that 
a_ j  =  a’  and  }  —  b *.  The  span-wise  velocity  at  the  surface  must  then  be 

v(x,X,y,t)  =  YljLi  ~r[ajy(X,y)  +  0(«)]c‘(*J'r-a’J<>  +  c.c. 

(2.66) 

+  Zj°=i  -^IM^y)  +  0(«)]e,(-k'I^t)  +  c.c. 
in  order  to  satisfy  Equation  (2.48).  Since,  to  lowest  order.  u0t  +  Vox  =  0,  an  expression  for 
the  surface  amplitude  is  readily  available:  the  replacement  of  the  lowest  order  velocity 
into  the  momentum  equation  yields 

OO 

1  K3 

oo 

+  E  T 9)  +  V)}ei{'k>x-*>t'  +  C.c.  (2.67) 

j= i 

A  solution  of  the  form  given  by  Equations  (2.65),  (2.66)  and  (2.67)  is  valid,  provided 
that  the  following  relation  holds  between  the  frequency  and  the  wavenumber: 

=  0.  (2-68) 

HJ’I 

which  gives  the  dispersion  relation  for  the  jth  mode,  the  positive  root  kj  corresponding 
to  the  shoreward-directed  wave,  and  the  negative  to  the  seaward  wave. 

The  solution  must  also  satisfy  a  compatibility  cond:'ion.  Since  the  linear  operator  C 
in  Equation  (2.60)  appears  in  every  order,  and  terms  of  lower  order  appear  in  the  inho¬ 
mogeneous  part,  secular  tjrms  arise.  It  is  an  artifice  of  having  truncated  the  expansion. 
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and  is  typified  by  the  possibility  of  blow-up  due  to  resonance.  This  resonance  condition 
for  jth  interacting  waves  is, 


kj  ±  •  •  •  ±  k2  ±  ki  =  0 

j  +  •  •  •  4-  —  0. 


(2.69) 


where  the  wavenumbers  and  corresponding  frequencies  obey  the  dispersion  relation  given 
by  Equation  (2.46).  In  the  scaling  adopted  in  this  study  the  0(kj )  =  0{u>j). 

The  resonance  condition  is  possible  for  certain  wave  systems.  For  isotropic  waves,  i.e. 
u>  =  u(k),  the  resonance  condition,  with  w(0)  =  0,  can  always  be  satisfied  if  uj'(k)  >  0, 
and  u"(k)  >  0.  These  conditions  cannot  be  satisfied  if  o/(k)  >  0,  u ;"(k)  <  0. 

To  prove  this,  we  note  that  the  dispersion  relation  u(k)  is  convex  downwards  if 
uj"(k)  >  0.  Hence, 


w(|k2  -  Ki  |)  <  w(kx)  +  u;(k2)  <  u>(k2  -f  Ki).  (2.70) 

Let  K3  =  K2  +  K1,  so  that  when  the  angle  between  K2  and  Kx  changes  between  0  and  tt, 
one  has  that  |«2  —  kx  |  <  K3  <  |k2  +  |.  Hence  the  frequency  W3  =  U3(k3)  will  change 

continuously  between  w(| k2  -  «x  j)  and  w(k2  +  kx),  and  will  coincide  with  u(ni)  +  u(k2) 
at  some  point.  H u/'(k)  <  0,  however,  we  have  u>(kx)  +  oj(k2)  >  u>(k2  +  «x)  >  u(n3)  and 
at  such  point  coincidence  is  impossible. 

Since  the  dispersion  relation  for  gravity  water  waves  is  such  that  u'(n)  >  0,  and 
u"{k)  <  0,  perfect  coincidence  is  not  possible.  At  most  we  expect  what  we  shall  refer  to 
as  “weak  resonance”.  Furthermore,  we  shall  restrict  our  attention  to  the  special  weakly- 
resonant  triad  case  in  which  k2  =  2k\  -  6,  u2  =  2wx,  where  the  detuning  parameter 
f)  <  0.  The  compatibility  condition  is  that  solutions  of  the  a,+1  be  orthogonal  to  the 
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solutions  of  order  a1  and  below,  so  that  resonant  solutions  are  avoided: 


jk i  fXo+2v^  +t]ki 


3*i  ( 
9*  Jx 


e±,3klI(G:  +G*)dx  =  0,  where;  =  1,2,3..., 


(2.71) 


starred  quantities  conjugated.  The  case  of  quartet  interactions,  to  lowest  order,  appears 
in  Appendix  B,  and  is  a  straight-forward  extention  of  the  ideas  presented  here. 

The  justification  for  using  a  small  number  of  modes  comes  from  field  data.  Figure 
2.3  suggests  that  most  of  the  energy  in  the  waves  is  found  in  the  first  few  harmonics. 
This  situation  is  quite  typical.  The  figure  also  shows  the  shifting  of  energy  from  lower 
frequencies  to  higher  ones  as  the  wave  travels  shoreward  over  a  decreasing  water  column 
depth. 


Figure  2.3:  Energy  for  shallow  water  waves  in  the  Southern  Baltic  Sea:  h0  =  6.0m 
- ,  ho  =  2.0 m - .  From  Druet  et  al.  [56]. 

Application  of  the  compatibility  condition  to  the  lowest  order  terms  in  Equation 
(2.63),  yields,  after  some  algebra,  the  evolution  equations  for  the  modes  in  Equation 
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(2.65),  and  Equation  (2.67) 

aix  +  iefDiEidi  -  iaFidiyy  +  iaDiSie~'^xa\a2  =  0 

a,2x  +  i£fD2E2a2  -  iaFiaiyy  +  ictD2S2e+'^xa1  =  0 

(2.  /  2) 

bix  -  iefD\E\b\  -f  iaF\biyy  -  iaD\S\e+'^xb\b2  =  0 
b2x  ~  ie fD2E2b2  +  iaF2b2yy  -  =  0. 

to  0(6/ X),  after  substituting  back  X  =  ax.  The  constants  are 

Dj  =  [2(1  -/J2^)]"1 

Ej  =  kj(  l-§/32^2) 

Fj  =  1/2  kj  (2.73) 

51  =  *^{*2 -*i +»!(%■  +  %)} 

52  =  2Atj/w2  +  2ul 

Equation  (2.72),  along  with  appropriate  boundary  conditions  determines  in  an  ap¬ 

proximate  way  the  ocean  surface.  The  incident  and  reflected  waves  are  decoupled  owing 
to  the  assumptions  and  restrictions  on  the  spatial  variation  of  the  bottom  topography. 


If,  on  the  other  hand,  the  longshore  sand  ridges  being  considered  were 

h(x,X,  y,T)  =  l  +  ef(x,X,y,T)  (2.74) 

the  resulting  modal  equations,  to  lowest  order,  would  be 

alx  -  isfDiEihai  +  iefDiEiPxbi  -  iaFialyy  +  iaDiSie~'^xa\a2  =  0 

a2x  ~  t£/^2^272a2  +  i£fD2E2fi2^2^x  ~  iaF2a2yy  +  iaD2S2e+t^xa\  =  0 

bix  +  isfD\Ei‘y1bi  -  iefD\E\^a\  +  iaF\biyy  -  iaDiSie+,^r6}&2  =  0 

b2iJri£fD2E2b2-i£fD2E2e~'2^x^2a'iFioiF2b2yy-ioLD2S2£~^:,:b\  =  0, 

(2.75) 
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to  0{f>l X),  with 

iki  /'2ir/ji  i 

Ij  =2^  J0  (/«  +  '2ikjfx  -  k]f  )dx 

(fxx  +  2ikjfx  —  k?j  f)e~2'iklX dx 

4  =  ^  j^hk\fxx  +  2ik]fx-k]f)e+2i^xdx.  (2.76) 

The  most  striking  difference  between  the  way  Equations(2.72)  and  (2.75)  describe  the 
surface  is  that  in  the  former  case,  the  bars  do  not  act  as  scatterers,  and  all  energy  in  the 
reflected  component  is  put  into  it  through  the  boundary  conditions. 

2.3  Internal  Wave  Case 

Shallow  water  weakly  nonlinear  interfacial  waves  appear  as  highly  coherent  groups  having 
well  defined  wavelength  and  are  observed  propagating  shoreward  on  a  density  stratifica¬ 
tion,  such  as  the  picnocline.  Their  crests  are  generally  oriented  along  isobaths  [57,58]. 
Their  wavelengths  range  from  200  to  1600  meters,  depending  on  the  depth,  which  can 
be  considerably  larger  than  the  local  water  column  depth.  An  estimate  of  the  energy 
contained  in  the  larger  ones  is  in  the  order  of  0.1  MJ/m2.  They  have  been  seen  to  appear 
twice  a  day  in  some  areas,  coinciding  with  the  tidal  cycle,  and  originate  mostly  in  places 
where  there  are  sharp  changes  in  the  bottom  topography,  such  as  on  the  edge  of  the 
continental  shelf. 

2.3.1  Internal  Wave  Hydrodynamic  Problem 

In  this  section  the  Hamiltonian  formulation  to  the  two-fluid  internal  wave  problem  is 
developed,  relying  on  Bowman’s  work  [50].  Illustrated  in  Figure  2.4,  define  the  domain 
fli  ss  R2  x  [ -H ,  77],  and  Q2  =  R2  x  [*7i  D],  The  lower  layer  (1)  has  a  uniform  density  px, 
and  the  upper  layer  (2)  a  density  p2  <  px.  The  fluid  is  subjected  solely  to  gravitational 
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forcing.  The  velocity  field  is  now  given  in  each  layer  by  (u.  w)t,  where  the  subscript 
refers  to  layer  1  or  2.  The  interface  between  the  two  fluids  is  given  by  z  =  ri(r.t)  and 
the  bottom  by  2  =  -H{t,T). 

,  _  i-  r _ 


2  = 


Figure  2.4:  Side  view,  interned  wave  problem. 

Assume  the  fluid  is  incompressible  and  irrotational  in  each  layer.  In  terms  of  a  scalar 
potential,  the  velocity  is  given  by 

(u  ,«:),  =  (2.77) 

From  conservation  of  mass,  Equation  (2.4),  the  equations  of  motion  within  the  fluid 
are 

A3d>i  =  0,  in  fi,.  (2.78) 

At  the  interface,  the  pressure  is  continuous,  hence  the  dynamical  boundary  condition 
is 


4>i,t  ~  "2  lV3<i>«|2  -SPiV^  at  2 


(2.79) 
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The  bottom,  which  is  assumed  impermeable,  has  a  normal  velocity  that  agrees  with 
that  of  the  fluid.  Thus 

4>i<2  =  -VH-^01  at  z  =  -H.  (2.80) 

The  kinematic  condition  on  the  interface  is  again  D =  0,  or 

<t>i,z  =  %  +  at  2  =  T).  (2.81 ) 

Finally,  we  make  the  simplifying  assumption  that  the  normal  velocity  disappears  at  the 
constant  air-water  interface: 

<t>2iZ  =  0,  at  z  =  D.  (2.82) 

2.3.2  Hamiltonian  Formulation  of  the  Internal  Wave  Problem 

The  conjugate  variables  in  this  case  are  rj  and  U  =  p2Vd>2  ~  The  Hamiltonian 

system  takes  the  form 

_  ,SE. 

6F 

Ut  =  -V(~).  (2.83) 

By  virtue  of  the  boundary  conditions  at  the  air-water  interface,  the  results  from  the 
previous  section  shall  be  exploited  to  arrive  at  a  regularized  Boussinesq  equation  for  the 
internal  wave  case.  The  potential  energy  is  simply 

v=  J  d2r^g{pA  -  p2)T)2.  (2.84) 

RJ 

The  total  kinetic  energy  is  the  sum  of  contributions  from  both  layers,  thus 
v  D 

A  =  pj  y  d2r  J  —  | V30J \2 dz  +  p2  J  d?r  J  —  |V3d>2|2dz  =  A i  +  A 2.  (2.85) 

RJ  ~H  K2  “*» 


The  kinetic  energy  in  the  lower  layer  is  given  by  Equation  (2.31).  We  need  only  to  figure 


out  A'2.  The  boundary  condition  given  by  Equation  (2.82)  can  be  exploited  to  find  A'2 
as  a  surface  integral,  using  Green’s  theorem.  Assuming  the  gradients  of  the  potential 


tend  to  zero  as  |r|  —  oo. 


A'2  =  -  J  d2rp2V<hV$2. 


(2.86) 


Define  the  pseudo-differential  operator  G  =  -kcoth(HkD),  which  comes  from  sat¬ 
isfying  the  boundary  conditions  on  the  interface  and  on  the  ocean  surface.  Adding  the 
expressions  for  A'i  and  A'2,  using  the  definition  of  U,  and  the  operator  G,  the  total 
kinetic  energy  is 

K  =  i  f  d2r{—(H  +  r,) U2  +  ^p-U-GU  +  U)2}  +  0(a3^2),  (2.87) 

i  J  0 1  pi  pi 


or  rearranging. 


K  =  l  f  dM^[(^  +  »?)-:^)U2  +  ^U-MU  +  ^£(V7/.U)2}  +  0(QV), 

2  J  p1  p\U  p\  p\ 

R2 

(2.88) 

where  M  =  ^  +  G=  p-  kcoih[HkD). 

Depending  on  the  size  of  D/X,  there  are  three  physically  distinct  possibilities: 

•  if  D/X  <  1,  then  U-A/U  =  0(a2 D/ X2),  and  M  %  l/D  —  D ^  .  For  such  case,  the 
terms  jjU-U  and  U  MU  balance  if  a2D/ X2  ~  1.  We  obtain  a  Boussinesq  system. 

•  for  D/X  ~  1  then  U-A/U  =  0(a2/X).  In  such  case,  if  aX2/D  ~  1,  we  obtain  the 
Intermediate  Long- Wave  equation. 

•  Finally,  D/X  >  1,  so  that  U  MU  =  0(a2/ A),  then  M  «  |fc|.  If  aA  ~  1,  the  result 
is  the  Benjamin-Ono  equation. 


Note  that  this  last  case  corresponds  to  a  very  deep  upper  layer,  lying  over  a  thinner  lower 
layer,  which  we  do  not  consider  relevant  to  the  problem  in  this  study. 
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By  substituting  the  expressions  for  the  potential  and  kinetic  energy.  Equation  (2.84) 
and  Equation  (2.88),  into  Equation  (2.83),  the  general  dynamical  equation  for  the  dy¬ 
namics  of  the  internal  wave  field  is  obtained: 

nt  =  -V-{[^-(ff  +  ri)  -  Su}  -  §V.{//(^)2U  + 

(2.89) 

U t  =  +  (Pi  -  P2)av}- 

The  result  from  linear  theory  can  be  recovered  by  neglecting  second  and  higher  order 
terms  in  Equation  (2.89).  The  solutions  proportional  to  exp {i(kx  -  wt)}  satisfy 


Thus, 


“1  =  +  e^ri)tHU 

ujU  =  k(pl  -  p2)gg. 


(2.90) 


<2  s  p  =  9-l~x  P2'[l  -  ^Hkcoth(kD)}. 


The  relevant  case  in  this  study  is  the  first  one.  The  Boussinesq  equation  is  then 


Vt  =  -V  {T-(ff  +  ij)U}-  £*V{j?(VJ?)2U  +  Itf2I>VV-U} 

U‘  =  -V-{^U-U  +  (Pl-p2)gV}. 

2.3.3  Regularization  and  Scaling 


(2.92) 


The  lowest  order  relations 


*  =  U} 


(2.93) 


U,  =  -V{(pj  -  p2)gr)} 

are  used  in  a  manner  similar  to  the  surficial  case  to  modify  the  troublesome  parts  of  the 
dispersive  terms  to  get  the  model  for  the  hydrodynamics  relevant  in  this  study: 


(2.94) 


th  =  -v-{£(ff-M)u}  +  |^v.[v(/rifc)] 

U‘  =  -V{^u-U  +  (Pi-P2)S7?}. 

The  scaling  is  the  same  as  the  surficial  case,  except  that  U  has  units  of  momentum. 
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Let  7  =  ^Pi-pPd  be  the  Boussinesq  parameter,  and  the  typical  thickness  of  the  lower 
layer  be  ho.  The  scaling  adopted  here  is 

»-■>/<■  i-g  r-i-  <2-93» 

Equation  (2.94)  is,  in  scaled  variables, 

qt  +  V-{Jr(/i  +  ar?)U}-d/J2^V.[V(%)]  =  0 

Ui  +  V{iU-U+7q}  =  0. 

Additionally,  the  span-wise  variables  are  scaled  to  reflect  the  weak  longshore  dependence 

of  the  waves: 

y  —  ol/2y  y-u  <—  a~l^2y  u.  (2.97) 

2.3.4  Slightly  Resonant  Interacting  Triads 

Once  the  equations  are  non-dimensionalized.  the  procedure  used  to  arrive  at  the  lowest- 
order  theory  is  the  same  as  was  done  for  the  surface  case.  In  this  presentation  the 
reflected  wave  shall  be  eliminated  from  the  outset. 

The  uniform  expansion,  after  cross  differentiating  Equation  (2.94),  is 

a0  :  Ci} o  =  0 

a1:  Ct) i  =  Gi(vo.uo,v0,G;x,X,y,t)  (2.98) 

a2:  £772  ,T}\,ui,vi,G;x,X,y,t) 

where 

£  =  3tt  -  -fdxx  -  d^-P^’1"  ,  (2.99) 

The  inhomogeneous  terms  Q\  and  Q2  involve  terms  that  appear  in  the  left  hand  side  of 


equations  of  order  lower  than  where  they  show  up.  Consider  the  lowest  order,  in  which 


Chapter  3 

The  Mass  Transport  Problem 


The  drift  velocity  is  the  second-order  steady  state  flow  that  is  created  by  the  passage 
of  overlying  water  waves  in  the  sediment-laden  boundary  layer  that  hugs  the  bottom 
topography.  The  boundary  layer  is  assumed  to  have  a  characteristic  thickness  6bi  <  ho- 
The  sediment  in  the  boundary  layer  shall  be  assumed  to  move  from  place  to  place  at 
a  rate  equal  to  the  drift  velocity.  This  chapter  shall  present  the  lowest  order  theory, 
leaving  the  details  of  the  higher  order  theory  to  Appendix  A. 

In  order  to  subsequently  compute  the  drift  velocity,  it  shall  be  required  to  find  the 
fluid  velocity  immediately  outside  of  the  sediment- laden  boundary  layer.  From  Equation 
(2.29)  in  scaled  variables,  the  shoreward  velocity  is  explicitly 

Ub  =  x-U(r  ,-M) 

=  u(r ,t)  -  02{—h[(huxx(r,  t))  +  a{hvxy(r,t))}  +  ±h2(uxx(rj)  +  avxy(r,  <))}, 

(3-1) 

and  the  span-wise  velocity 
Vb  =  j)-U(r,  —h,  t) 

=  v(r,t)  -  P2{-h[(huxy(r,  t))  +  a(/inyy(r,0)]  +  \h2{uxy{  r,t)  +  avyy{r,t))} 

(3.2) 

in  the  neighborhood  of  the  boundary  layer.  Neglecting  the  reflected  component,  the 
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In  the  boundary  layer  the  transverse  momentum,  vertical  momentum,  and  the  continuity 
equations  are  respectively 

Ut  +  u-Vu  +  iuu2  =  -iVp  +  i/Au  +  i/u*2 

wt  +  u-Vti;  +  wwz  =  -  ^Pz  +  9  +  t'Ati;  4-  uwzz  (3-4) 

V-u-)-u)z  =  0 

where  u  is  the  eddy  viscosity.  Across  the  boundary  layer  the  flow  velocity  changes 
from  zero  at  the  bottom  boundary  to  some  finite  value  characteristic  of  the  exterior 
inviscid  fluid.  The  derivatives  with  respect  to  z  of  any  flow  quantity  are  thus,  in  general, 

much  greater  than  those  with  respect  to  x  or  y.  Hence,  within  the  boundary  layer, 

|Vu|  <  |ti*|,  |V2u|  <  |u«|,  etc.  In  view  of  this,  it  may  be  concluded  that  the 
transverse  momentum  in  Equation  (3.4)  is  well  approximated  by 


u<  +  u-Vu  -|-  wuz  =  — Vp  +  vuzz. 

P 


(3.5) 


The  velocity  w  must  also  be  small.  The  continuity  statement  in  Equation  (3.4)  suggests 
that  the  boundary  layer  and  w  are  of  equal  order  of  smallness.  Therefore,  none  of  the 
terms  on  the  left  hand  side  of  Equation  (3.5)  can  be  neglected.  On  the  other  side  of 
the  equation,  it  is  expected  that  i/u„  be  of  comparable  size  to  the  inertial  terms.  The 
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magnitude  of  the  inertial  terms  is  represented  by  the  size  of  u-Vu,  hence  the  balance  is 
such  that  0(u-Vu/i/u„)  =  1  when  the  Reynolds  number  is  sufficiently  large,  i.e.  so  it 
is  mostly  true  in  the  whole  boundary  layer.  If  y/gh0  is  representative  of  the  magnitude 
of  the  velocity  u  and  A  represents  a  distance  in  the  transverse  direction  over  which 
u  changes  appreciably,  then  {y/gho)2 /A  =  O(u-Vu).  Since  bbt  is  the  boundary  layer 
thickness,  i/y/gho/b2[  is  a  measure  of  vwZ2.  Thus, 

0{62blR/ A2)  =  1,  where  R  =  (3.6) 

The  dimensionless  constant  R  is  the  Reynolds  number.  Equation  (3.6)  implies  that 
6bi  ~  A R~1/2  as  R  — ►  oo,  which  suggests  that  the  boundary  layer  concept  improves  as 
R  — »  oo,  and  that  6bl  «  A ll2vll2.  In  this  study  A  is  large  but  of  finite  length.  It  is 
assumed  that  the  boundary  layer  does  not  change  significantly  as  a  function  of  wave 
frequency,  enabling  the  replacement  of  A  by  h0  in  R ,  so  that  R  =  y/gh0h0/i/,  arriving, 
then,  at  a  working  definition  for  the  boundary  layer  thickness 

hi  =  \Jv/ho(gho y/2,  (3.7) 

which  shall  be  non-dimensionalized  by  dividing  by  ho.  In  this  scaling,  it  is  implied  that 
the  size  of  the  Reynolds  number  and  the  boundary  layer  thickness  are  controlled  mostly 
by  the  viscous  effects,  i.e.  the  size  of  u. 

To  get  an  estimate  of  the  size  of  w,  the  continuity  condition  in  Equation  (3.4)  suggests 
that 

0(u>*/Vu)  =  6biho/X  =  3y/gh^R~1/2  (3.8) 

thus 

w  =  O(0y/^R-'12).  (3.9) 


With  Equation  (3.9)  in  hand,  it  can  be  inferred  from  the  vertical  momentum  balance 


0  I 


that 

g  +  —  =  0(Su).  (3.10) 

P 

hence  pz  =  O(Sbi),  i.e.  the  pressure  is  approximately  constant  throughout  the  layer. 

We  are  now  in  the  position  of  estimating  the  balance  of  terms  in  Equation  (3.4).  We 
approach  this  in  stages.  Firstly,  to  make  the  system  in  the  boundary  layer  consistent 
with  those  in  the  inviscid  fluid,  we  adopt  the  inviscid  scaling.  The  equations  are  now 

j3ut  +  q/3u-Vu  +  awuz  =  -§Vp  +  -^=-[/32Au  +  i/uzz] 

(3wt  +  a/9u-Vu)  +  awwz  =  -}rp2  +  \la-\ — ,  *  [/32Aw  +  vw,z 1  (3.11) 

9"  oho 

/3Vu  +  wz  =0. 
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A  locally  flat  bed  has  been  assumed.  In  contrast,  suppose  that  the  bed  had  some 
finite  curvature  A',  say.  This  would  change  the  vertical  momentum  balance  in  Equation 
(3.14)  to 

pn  =  A'0(u2).  (3.15) 

but  the  pressure  change  across  the  layer  is  still  of  0(f>bi ).  so  we  are  justified  in  the 
assumption  that  the  bed  be  locally  flat. 

The  following  boundary  data  is  used  to  solve  Equation  (3.14): 

u  =  v  =  w  —  0  at  n  =  0  (3.16) 


and 


u^Ub 

v  —  Vb,  n  — *  oo. 


(3.17) 


The  velocity  (Ub,  H)  immediately  outside  of  the  layer  gives  rise  to  the  following  pressure 
gradients: 


(3.18) 


-qPx  =  PUbt  +  (*t3(UbUbl  +  aVbUby) 

i py  =  PVbt  +  a3(UbVbx  +  aVbVhy). 

We  have,  thus,  all  the  required  information  to  solve  for  the  velocities  in  the  boundary 
layer.  Performing  the  usual  expansion 


u  =  Uq  +  otu i  ■  •  • 


v 


v0  +  avi  ■  ■  ■ 


(3.19) 
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the  lowest  order  equations  are 

(3uot  —  Uonn  —  fiUobt 

dv0t  ~  UOnn  :  r^^Obt 

(3.20) 

POn  =  0 

fiuox  +  Won  =  0. 

The  following  order  is  given  by 

@U\t  —  U\nn  =  —(3UoUox  —  WqUOu  +  PUobUox  +  Ult 

0Vlt-Vinn  =  ~/3u0V0x  -  W0V0n  +  PUobVox  +  dVibt  (3.21) 

f3uix  +  u)in  =  -f3v  0y, 

and  shall  be  addressed  in  Appendix  A.  Note  that  terms  such  as  [fu,  etc.,  have  been 
dropped.  The  goal  is  to  compute  the  drift  velocity  and  these  terms  do  not  contribute  to 
the  steady  part  of  the  drift  velocity  since  the  external  flow  is  time  harmonic. 

A  solution  of  Equation  (3.20)  of  the  form 

2 

ut  =  ^2  a,Pi(xi  V'  n)el(k]X~U]t^  +  c.c.  (3.22) 

j= i 

subject  to  the  boundary  conditions  given  by  Equations(3.16)  and  (3.17),  is  found  by 
integrating  Equation  (3.20).  The  same  procedure  is  used  to  obtain  v.  The  result  is 

iio  =  Cja;(l  -  +  c.c. 

vo  =  *Ei=i^(^fcVi/2“Q°j»/*i)(l-e_BA,)e,'(,i'<-w^)  +  c.c.  (3-23) 
Wo  =  ~  *Aj  -  e-n^>)/AJe,(fcJx_ubt)  +  c.c. 

where  Aj  =  (1  -  2.  The  vertical  velocity  w  was  found  by  integrating  the  conti¬ 


nuity  equation. 
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3.2  The  Drift  Velocity 

In  this  section  we  follow  closely  Longuett-Higgins’  study  on  mass  transport  by  oscillatory 
flows  [25].  Define  the  time  average  of  the  quantity  A  as 

U)  ft+  uT  1  ft+T 

{A)  =  —  /  A{s)ds  =  -  A(s)ds.  (3.24) 

2tt  Jt  Tjt 

The  drift  velocity  shall  be  the  time  average  displacement  rate  of  a  fluid  particle. 

As  Stokes  noted  [26],  the  drift  velocity  is  second  order  in  nature.  That  is,  if  u(r,  z,t)  is 
the  Eulerian  velocity,  and  the  motion  is  periodic, 


u(r,z,t  +  r)  =  u(r,z,t), 
and  expressible  as  the  asymptotic  series 

XL  —  U0  +  G^tti  +  a2«2  +  •  '  • 


(3.25) 


(3.26) 


then  (iio)  =  0  is  a  statement  of  this  assumption,  i.e.  the  lowest  order  steady  state  current 
is  zero.  Let  U(r0,2o,io)  denote  the  Lagrangian  velocity,  or  velocity  of  a  fluid  particle 
at  t  =  t0  with  position  (r0, 2o)-  Then,  the  displacement  of  the  particle  from  its  original 
position  to  some  other  position 


(r,*)  =  (r0,20)+  f  U(r0,zo,i)dt 
Jto 


(3.27; 


at  time  t.  It  follows  that 


U(r,z,<)  =  u[(r0,2o)  +  /  U{r0,zo,i)di,t} 

Jt0 


(3.28) 


which  can  be  formally  expanded  in  a  Taylor  series 


U(r,  z,t) 


u(r0,  *o,  t)  +  ft[  Ud<'V(ro-2o)^(rOi  t) 

+  |{//0  Udt}T.H{To,Zo)u(r0,z0,t)-f;o  Udi  +••• . 


(3.29) 


Here  H  stands  for  the  Hessian,  and  both  the  Hessian  and  gradients  have  subscripts  to 
remind  us  as  to  where  they  are  to  be  evaluated.  The  steady  state  Lagrangian  velocity 
is  in  fact  akin  to  the  drift  velocity.  An  approximate  expression  for  it.  in  terms  of  the 
Eulerian  velocity  in  the  boundary  layer,  can  be  obtained  by  expanding 

U  =  Uo  +  a'Ui  +  a2U2  +  ---  (3.30) 


and  substituting  Equations  (3.30)  and  (3.26)  into  Equation  (3.29),  order  by  order.  After 
time-averaging  we  obtain 

O(o°)  :  <Uo>=  (u0)  =  0 

Ofa1):  (Ut)=  (Ul)  +  </‘  Uodf-Vuo) 

(3.31) 

0(a2)  :  <U2)  =  (u2)  +  (f  Uidt-Vuo)  +  (f*  U0dt-VUl)+ 
i({/‘  u0di}T  Hu0-f‘o  u0 dt). 

The  drift  velocity  is  then 

(U,V)  =  a1Ul  +a2U2  +  --.  (3.32) 


We  shall  compute  the  drift  velocity  to  lowest  order  in  this  chapter,  and  leave  consid¬ 
eration  of  the  higher  order  contribution,  which  can  be  rewritten  as 


(U2)  =  (u2)  +  (f‘  Uidt-Vuo)  +  (/'  f'  uodt'-V  u0df- V  u0)  +  (fl  u0dt-Vuj)+ 

i({/‘  uodt}T  H u0-  //o  u 0dt), 


(3.33) 


to  Appendix  A.  Expressing  the  0(a)  in  component  form,  after  weak  y  dependence  scaling 
has  been  adopted,  the  drift  velocity  is 


—  (til)  +  (ft  UodtUox)  +  (/*  WodtUQn) 
v  =  (vi)  +  (f*  U0dtj0x)  +  (/‘  w0divQn). 


(3.34) 
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3.3  The  Mass  Transport  Equation 

3.3.1  Remarks 

In  regions  where  sand  ridges  are  found  the  mean  slope  of  the  bottom  is  very  slight. 
Hence,  downslope  gravitational  transport,  important  in  the  shoaling  region,  plays  an 
undiscernible  role  in  shaping  sand  ridges  in  the  deeper  end  of  the  shelf.  The  type  of 
sediment  transport  we  have  in  mind  is  primarily  suspended  load.  The  sediment  drifts  as 
a  result  of  the  unclosed  orbital  paths,  resulting  from  the  asymmetry  in  the  nonlinear  wave 
motion.  In  what  follows  the  fluid  wave  field  shall  be  assumed  to  be  entirely  represented 
by  the  incident  wave.  Further,  it  is  assumed  that  the  viscous  boundary  layer  is  sediment¬ 
laden,  composed  of  cohesionless,  rarely  interacting,  sand  particles. 

The  sediment  concentration  p,  in  coastal  environments  has  a  very  weak  influence  on 
the  fluid  flow  [59].  Typical  values  for  the  concentration  are  p  ~  10-3  -  10_4ppm,  and 
it  shall  be  assumed  that  this  is  the  situation  thoughout  the  shelf.  Chapalain  [59]  and 
Bokzar-Karakiewikz  et  al.  [60]  concluded  that  time  independent  and  vertically  uniform 
parameters  of  eddy  viscosity  and  eddy  diffusivity  are  adequate  in  providing  satisfactory 
accuracy  for  sediment  morphology  models  on  the  shelf.  In  this  study  we  shall  adopt  a 
very  simple  model,  found  in  [3],  for  the  sediment  concentration. 

The  equation  of  continuity  for  the  sediment  concentration  is  the  advection-diffusion 
equation 

p(  +  V-(up)  +  [(w  -  vf)p\n  =  0,  (3.35) 

where  vj  is  the  sediment  “fall  velocity" and  n  =  {z  +  h(X,y,T))/6bi-  Assume  that,  appart 
from  random  fluctuations,  u  and  p  do  not  vary  much  over  small  transverse  spatial  scales 
so  that  the  second  term  of  the  above  equation  may  be  neglected.  In  light  of  this,  the 
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sediment  concentration  changes  at  a  rate  dp/dn  proportional  to  the  vertical  flux.  Hence. 

wp=-lPn  (3.36) 

where  7  is  the  diffusivity  constant. 

The  flux,  which  is  the  product  of  the  concentration  and  the  velocity,  can  be  split  into 
a  time  dependent  part  Cl  and  a  time  independent  part  Cm.  Boczar-Karakiewicz  et  al. 
[9]  found  that  in  the  sand  ridge  areas,  the  ratio  C(/Cm  =  O(10~2)  for  the  off-shore  case. 
This  situation  shall  be  assumed  to  apply  throughout  the  shelf,  so  that  the  sediment 
concentration  shall  be  represented  solely  by  the  time  independent  part  in  this  study. 
Employing  this  assumption  and  substituting  Equation  (3.36)  into  Equation  (3.35),  the 
equation  for  sediment  concentration  is  now 

lPn  +  vJP-Q-  (3.37) 

The  boundary  condition  may  be  taken  as 

JJK  =  P(r)’  ,W8) 

where  P( r)  is  akin  to  Svendsen’s  [61]  “pick-up  function’’,  which  incorporates  such  effects 
as  the  degree  of  wave  asymmetry  and  skewness  of  sediment  flux,  and  a  Shield’s  parameter, 
which  sets  a  threshold  fluid  velocity  at  which  sediment  will  be  picked  up,  based  upon  the 
sediment  particles’  buoyancy  and  geometry,  and  the  fluid’s  velocity  field  and  viscocity. 
The  pick-up  function  P(t)  is  obvoiusly  an  empirically-derived  function. 

Solving  Equation  (3.37),  the  sediment  concentration  takes  the  form 

p=P(  r)e~vn,  (3.39) 

where  a  —  Vf/ 7.  The  fall  velocity  v/  is  species-dependent.  It  is  either  measured  or 
estimated  by  calculating  the  balance  of  drag  to  buoyant  forces  for  a  particle  falling 
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freely  into  a  static  fluid.  The  diffusivitv  constant  7  is  hard  to  estimate.  Sedimentologists 
usually  measure  its  value  in  the  field.  The  pick-up  function  P( r)  is  an  empirically-derived 
function. 

3.3.2  The  Transport  Equation 

For  the  sake  of  clarity,  the  mass  transport  equation  shall  be  derived  assuming  transverse 
dependence  in  the  x  direction  only.  The  generalization  to  variations  in  y  will  follow  in  a 
straight-forward  manner. 

Let  T  €  [0,oc)  and  €  721  x  [/i(T),£],  where  £  >  h(T)  +  6bi ,  be  the  boundary  layer 
time-space  domain,  and  consider  a  differential  “volume”  element  in  such  a  domain,  as 
shown  in  Figure  3.1,  which  on  the  bottom  is  bounded  by  the  ocean  topography,  and  on 
the  top  by  a  flat  lid  z'  =  (.  By  examining  the  balance  of  mass  in  this  differential  volume 
a  transport  equation  can  be  developed. 

It  is  assumed  that  the  sediment  concentration  p  is  entirely  negligible  for  z'  >  (,  and 
moves  on  fast-time  scales.  In  what  follows  p  :  fir  ' — *  721.  The  sediment  concentration 
and  drift  velocity  are  thought  to  be  C'ffix).  and  the  bottom  topography  h  G  C'(T), 
and  piece-wise  linear  in  fij. 

The  mass  flux  per  unit  length  at  x  in  a  time  interval  [T,T  +  AT]  is  given  by 
/•r+Ar  k  rT+&T  k 

/  dr  dz'p(x,z')U{x,z')  =  /  dr  dz'M{x,z').  (3.40) 

Jt  Jt  A(x.r) 

Consider  a  portion  of  the  region,  say  [x,x  +  Ai]  in  a  time  interval  [T,  T+AT].  Since  mass 

cannot  spontaneously  vanish  or  be  created,  the  net  amount  of  sediment  between  point  x 

and  x  -I-  Ax  must  be  compensated  by  a  change  in  the  concentration  of  the  sediment,  or 

by  a  topographical  change  in  the  bottom  surface.  The  flux  difference  in  the  space  and 
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A T  go  to  zero,  we  have,  on  the  right  hand  side. 


Un^T-oA,-o  shf  //+ A'  *  JS+AT>  d -’V 

,  ,\  A(5.T)+At^A> 

^mAr-.oAx-.o  AxA T  ^  fh(i,T)  T  dz'p(t,:')  = 


p(x,h(x,T))^fl, 


and  on  the  other  side  of  the  equation. 


tomAr— oAr—o  Ax  At  It+~^T  dT{f^(  x+b,x,T)dz'  M(x  +  Aar'z')  ~  fh(x.T)dz'AI(x-  -')}  = 


limAr— oAx— o  AxAt  { J^(l|Ai,r)  dz/[A/( x,  z;)  4-  Ai^(i,:')  +  •••]- 

4.d  =  ^mAr—oAx— o  SZXrAx4C(l>r)+Ai^  dz'^(x’2,)  = 

()x 


fh(x,T)  dz'^fc(X'  Z>) 


Therefore,  the  mass  transport  equation  is 

d/i(x.T)  A"  d 


dT  p(x 


,  h(*, T))0x  /Mx.T)  **’ ^  (3’47) 


where  A''  is  a  constant  of  proportionality.  Since  the  boundary  layer  is  assumed  very  thin, 
we  may  define  the  mass  transport  flux  as 


p  =  fSb‘  p(x,z')U{x,z'))dz' 
Jo 


fVbl 

V  =  /  p(x,z')V{x,z'))dz' 

Jo 


(3.48) 


so  that  the  transport  equation  now  reads 


dh(x,T )  A' 
dT  ~  p /*■ 


(3.49) 


The  generalization  of  Equation  (3.49)  to  one  more  space  dimension  is 


dh(x.y,T)  _  A' 

dT  ~  pQifiz  +  l/y)' 


(3.50) 
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where  p  and  v  are  the  shoreward  mass  flux  and  the  longshore  mass  flux,  respectively. 
Note  that  when  weak  y  dependence  scaling  is  adopted  in  Equation  (3.50).  the  longshore 
mass  flux  shall  be  0(a)  smaller  than  all  other  quantities  in  the  equation. 

In  the  remainder  of  this  study,  it  shall  be  assumed,  for  simplicity,  that  the  sediment 
concentration  is  constant  and  equal  to  p0  in  the  boundary  layer.  In  terms  of  Equation 
(3.34),  and  upon  use  of  Equation  (3.48),  the  calculation  of  the  mass  flux  components,  to 


lowest  order,  are 


where 


f*  -  — ~zr — +  L  — 3 — J2j  +  c.c. 


JTi  *iai 


(3.51) 


+e~^u [cos Ojtu  -  sin  -  d<Tj{a36bl  +  1)]  (3.52) 


=  ^(1/2 -^w)  +  e-M«/4 

-e~<T}^bt[l  4-  6bl<7j]  cos (TjSbi  4-  2e~<7j^bl  sin  a: 


(3.53) 


for  the  shorew'ard  mass  flux,  and 


JL,  iC2a*a,v  , 

^  =  E  -- J  Jj  +  O(03)  +  c.c. 


j=i  uiai 


for  the  longshore  directed  mass  flux,  with 


(3.54) 


Jj  =  <TjSu-  1  — -(1-e  )  4-  e  <t>^w(cos<7j6(,j  -  sin  (Tjfbl) 

4-/3Aj[|(  1  +e-2°>6«)  +  e-X’Sb‘(t6bloJ/2-  1)]. 


(3.55) 


The  quantities  In,  I2j  and  J\  are  plotted  parametrically  in  Figures  3.2,  3.3.  and  3.4. 


Figure  3.2:  Variation  of  2n,  with  6^  =  1.0  fixed. 


Before  proceeding,  two  important  remarks  are  in  order.  Firstly,  it  is  noted  that  the 
bottom  needn’t  be  slightly  perturbed  to  initiate  the  development  of  bars.  Even  a 
flat  bottom  will  eventually  develop  bars  given  that  all  the  conditions  are  right.  Sec¬ 
ondly,  we  are  now  in  the  position  to  justify  the  two-time  scale  solution  of  the  sur¬ 
face/bottom  system.  The  ratio  of  the  magnitude  of  the  time  rate  of  change  of  the 
bottom  to  the  Eulerian  velocity  leads  in  a  straight-forward  manner  to  the  conclusion 
that  t/T  =  0(a)0(f>bi)0(p)  ~  O(10-7),  assuming  that  the  boundary  layer  thickness  is 
typically  O(10-2/io)  and  the  sediment  concentration  is  O(10-4)  ppm. 


Chapter  4 

The  Complete  Model: 
Mathematical  Analysis 


After  summarizing  the  model,  this  chapter  shall  be  devoted  to  the  details  of  formal 
and  analytical  results,  primarily  relevant  to  discerning  the  behavior  and  structure  of  the 
surface  system.  Our  concentration  on  the  surface  system  is  motivated  by  the  fact  that 
much  is  known  already  about  the  type  of  equation  represented  by  the  mass  transport 
equation,  whereas  the  surface  system,  insofar  as  we  can  tell,  is  a  new  mathematical 
equation.  Many  of  the  surface  system  results  presented  here  are  actually  not  applicable 
to  the  oceanic  setting  in  which  this  equation  was  derived;  additionally,  some  the  results 
are  not  entirely  new.  However,  since  the  surface  system  is  interesting  in  its  own  right, 
and  it  could  serve  as  a  model  for  other  physical  processes,  the  results  are  still  of  value. 

For  the  surface,  after  replacing  the  multiple  scale  expansion  X  =  ax , 
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-  iK\aXyy  +  il\3f(x,  y)ax  +  iK5e~^xa\a2  =  0 

a2x  -  iK2a2yy ih'4f(x,y)a2  +  iK6e+l^xa1  =  0 
6lr  +  iLibXyy  -  iL3f{x,y)bi  -  iL5e+t^Ibi{b2  =  0 

b2x  +  iL2b2yy  -  iL4f(x,y)b2  -  iL6e~li>Ib\  =  0 

(4.1) 

ai(i  =  0, y)  =  A\{y) 

a2(x  =  0,  y)  =  A2(y) 

bx(x  =  M,  y)  =  Bx(y) 

b2(x  —  M,y)  =  B2{y) 

plus  appropriate  boundary  conditions  on  y  =  0  and  y  —  N .  The  K  and  L  coefficients 
are  0(a,e),  and  are  implicitly  given  by  Equation  (2.72)  and  Equation  (2.73).  If  the 
boundary  conditions  B\  and  B2  of  the  reflected  wave  are  small,  the  reflected  component 
is  negligible.  Assuming  this  is  the  case,  the  surface  system  is  then 

aXx  -  iKxaiyy iK3f{x,y)ai  +  iKse~'^xa4a2  =  0 
a2x  -  iK2a2yy  +  iK4f(x,  y)a2  +  i A'6e+,<5xa?  =  0 

(4.2) 

ai(x  =  0,y)  =  .Ai(y) 

a2(x  =  0,y)  =  A2(y). 

Although  the  linear  part  of  the  surface  system  is  identical  to  its  counterpart  in  the 
Nonlinear  Schrodinger  Equation,  the  nonlinear  terms  endow  the  surface  system  with 
properties  and  behavior  much  unlike  the  Schrodinger  Equation.  The  bottom  evolution, 
is  given  by  Equation  (3.50): 

&h<-,y,T)  =  %(nx  +  vy)  (43) 

h(x,y,  0)  =  H(x,y). 

Equation  (4.2)  and  Equation  (4.3)  comprise  the  full  model.  In  Chapter  3  we  gave 
an  estimate  of  the  time  discrepancy  for  the  evolution  of  the  surface  and  bottom.  This 
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discrepancy  suggests  an  effective  decoupling  of  the  fluid  and  sediment  problem,  which 
permits  an  iterative  solution  to  the  full  model.  Specifically,  begining  with  some  initial 
bottom  configuration  H{x,y ),  the  field  on  the  surface  is  solved;  the  flux  velocity  J 

is  calculated,  and  the  bottom  is  updated  using  the  mass  transport  equation.  With  this 
new  bottom,  the  fluid  quantities  are  solved  for  again  and  the  process  is  repeated  until 
some  T  final. 

The  conditions  for  which  the  surface  system  and  the  mass  transport  equation  are 
stable  enable  us  to  discern  the  conditions  for  the  overall  stability  of  the  iterative  solution 
of  the  model.  At  this  point  in  time  the  evidence  of  the  surface  system's  solution  stability 
comes  from  numerical  calculations.  Until  the  issue  of  stability  of  the  surface  system  is 
studied  in  detail,  it  shall  be  assumed  the  solutions  are  stable  and  in  what  follows,  and 
proceed  to  find  conditions  to  be  met  by  the  mass  transport  equation  so  that  the  overall 
iterative  procedure  is  stable. 

The  mass  transport  equation,  from  Equation  (3.50),  is  of  the  form 
dh  x^y,  T)  _  Vy)  jn  7£2  with  T  G  [0,  oe ) 

h(x,y,  0)  =  : H(x,y ).  (4.4) 

The  properties  of  this  quasi-linear  hyperbolic  equation  are  well  known  [62],  and  the  ex¬ 
istence  and  uniqueness  of  solutions  is  well  established:  provided  the  initial  condition 
h(x,y,  0)  =  7i(x,y)  is  at  least  in  the  C1  class  of  functions,  and  the  characteristics  are 
nowhere  parallel  to  the  manifold  on  which  the  initial  data  is  prescribed,  we  have  either 
solutions  in  the  weak  or  strong  sense,  i.e.,  smooth,  or  shock-like.  A  shock  solution  can 
either  be  prescribed  as  initial  data,  or  can  occur  at  some  later  time  when  the  character¬ 
istics  cross  in  space-time.  Regarding  the  problem  in  sedimentary  transport,  a  shock-like 
solution  would  make  little  sense  as  a  solution.  Since  the  possibility  of  such  an  outcome 
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exists,  it  is  worth  reviewing  the  conditions  for  the  formation  of  a  shock,  and  its  relevance 
to  the  problem  of  sedimentary  structure  formation.  Equation  (4.4)  may  be  recast  as 


dh(x,y,T) 

dT 

h(x,y, 0) 


Kidhhl+dhhy) 


H(a :,y) 


(4.5) 


assuming  that  the  indicated  differentiations  can  be  performed.  If  in  the  above  system  we 
identify  c  =  {-jjfc-sjfa)  =  (<~i > c2 )  as  propagation  speeds,  we  may  reinterpret  the  problem 
in  terms  of  simple  wave  dynamics.  Assuming  the  solution  is  wave-like,  it  may  be  inferred 


h  -  H( x  -  Tc\(h))  -  0(a) 


since  i/y  =  0(a)  the  second  term  on  the  right  hand  side  of  Equation  (4.5)  affects  the 
outcome  very  minimally.  Assume  that  H  is  differentiable.  Using  the  implicit  function 
theorem 


hr 

hx 


V! 

1  +  H'clhT 

W 

1  +  H'cxhT' 


(4.7) 


It  is  evident  from  this  pair  of  equations  that,  for  c\h  >  0,  if  H'  >  0  for  all  x,  both  hj  and 
hx  shall  remain  bounded  for  all  time.  On  the  other  hand,  if  cj*  <  0  at  some  point,  hj 
and  hx  shall  diverge  as  1  +  H,Ci/lT  —  0.  The  situation  is  the  reverse  if  c\h  <  0,  of  course. 
In  the  sediment  problem,  as  is  evident  from  Figures  4.1  and  4.2,  it  is  typical  for  Cj  to  be 
oscillatory  in  nature.  Hence  the  characteristics  have  wave-like  spatial  dependence. 

We  need  to  reiterate  the  issue  at  hand:  as  the  waves  shoal,  they  “see”  a  bottom  that 
is  essentially  fixed.  Only  after  the  passage  of  many  waves  do  we  expect  to  see  changes  on 
the  bottom  topography.  Insofar  as  the  solution  of  the  model's  system,  we  are  concerned 
with  the  issue  of  stability  in  the  iterative  solution  of  the  surface/bottom  system,  which 
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Figure  4.1:  ci(T  =  0)  for  bottom  f(T  =  0)  =  O.Olz.  0  =  0.0*,  s  =  0.2,  a  =  0.1, 

ta>i  =  1.2. 

comprises  the  full  model.  We  think  of  H  as  an  entirely  new  initial  condition  as  input  to 
the  conservation  law  at  each  value  of  T.  We  may  ask  then,  when  do  the  characteristics 
cross  so  that  the  solution  is  no  longer  valid?  In  what  follows  we  think  of  H  as  an  entirely 
new  initial  condition  and  T  as  the  time  between  any  two  iterates  of  the  full  model.  Set 

1  +  WWT  =  0.  (4.8) 


If  there  is  crossing  of  the  characteristics,  it  shall  occur  at 


T  = 


1 


H'clk 


(4.9) 


Figure  4.2:  C\ (T  =  0),  when  f(T  =  0)  =  O.Oly,  0  =  0.08,  e  =  0.2,  a  =  0.1,  u;j  =  1.2. 


For  the  two-dimensional  case 


c\h  =  -7iai/i  W  _  l?a\h{h), 

1  j  = 

/,  =  i-§, nkjh2 

thus,  crossing  occurs  when 

T=  — — L _ . 

'X'hialfiW  +  ')2alf2(h)\ 


(4.10) 


(4.11) 


By  assumption,  H1  =  O(e)  =  O(a).  Since  are  all  0(1),  and  fij  =  O{01^2) 
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then  7j  =  0(j33^2).  Thus,  an  estimate  for  the  time  at  which  crossing  may  occur  is 

H'12) 

which  can  be  quite  a  large  interval  assuming  that  |a,|  remains  bounded  and  less  than 
one.  This  estimate  applies  to  the  three-dimensional  case  reasonably  well  since  y  variation 
is  0(a)  smaller  than  x  variations.  Equation  (4.12)  gives  an  upper  bound  on  the  time 
intervals  between  each  iteration.  Recall  that  we  are  considering  each  iterate  as  a  new 
initial  condition,  and  that  the  drift  velocities  will  be  different  at  each  time  step.  Hence, 
if  the  upperbound  in  T  given  by  Equation  (4.12)  is  never  reached  -i.e.,  the  model's 
assumptions  are  never  violated  -  the  iteration  procedure  yields  a  stable  result. 

4.1  The  Surface  Equations 


As  mentioned  earlier,  the  surface  system  is  a  new  equation  whose  properties  and  structure 
are  presently  unknown.  In  what  follows,  we  present  what  we  have  been  able  to  discern 
thus  far  about  the  structure  of  this  new  system. 

4.1.1  Hamiltonian  Structure 


The  following  scaling, 

i  -  /f5(A'e£0)l/2i  y-K'AK tE0)‘l'y  A Q  -  Ki(K< 

14. IJ) 

V2d>  Jlib 

U  (A'.Eo)1'2  V  (KsEo)x!2  ’ 

shall  be  adopted  in  what  follows  in  order  to  facilitate  the  derivation  of  the  surface  equation 
Hamiltonian  structure.  For  a  flat  bottom,  the  system  is  thus 


ux  -  iK\uyy  +  ie  x^Iu*v  =  0 
vx-iK2vyy  +  ie+l*Qlu2/2  =  0 


(4.14) 


Assume  compact  support  in  y  for  the  dynamical  variables.  Define  the  Lagrangian  density. 

£  =  iuxu*  -  iu*u  +  ivxv*  -  ivxv  -  3f?[(u')2re~‘A^x]  _  A'i|u2|  -  A'2|c2|  (4.15) 

where  stands  for  “the  real  part  of”,  and  the  canonical  momenta 


nx  = 

dc  _ 

duz 

iu * 

n2  = 

dc  _ 

dvx 

iv* 

HI  = 

dc  _ 

dui 

iu 

ll 

G 

dc  _ 

(Jv* 

iv. 

(4.16) 


The  requirement  that  L  =  /  Cdy  be  stationary,  yields  the  Euler- Lagrange  equations, 
which  in  turn  lead  to  (4.14): 

=  '“*  +  =  0  ^  ^ 

=  i»r  +  A>„  -  e+iA«*«!/2  =  0 

and  its  complex  conjugates.  The  Hamiltonian  H  and  its  density  H  are  given  by 

H  =  J  Hdy 

H  =  ft[(u*)2re-,A<3l]  +  h\\uy\2  +  A'2|t’y|2.  (4.18) 


(4.18) 


Note  that  the  Hamiltonian  is  not  conserved,  i.e.. 


dL  dH  ,  n 
dx  ~  dx  *  °’ 


(4.19) 


except  when  t\Q  =  0.  The  Hamiltonian,  in  terms  of  the  conjugate  momenta  is 

n=  -i»[n?n;e-iA<?i]  +  A',|nly|2  +  a'2|  n2y  |2.  (4.20 ) 

The  Hamiltonian  admits  a  Poisson  structure:  Defining  the  Poisson  bracket  as 


r.  [JtdAdB  dA  OB  f  ,  .dA  dB  dA  dB . 

I-4- B>  55  J ^anT  -  anTa^1  +  cc' + J  ~  + 


dA  dB  dA  dB , 


c.c.  (4.21 : 
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so  that  the  evolution  of  a  dynamical  variable  .4  is  given  by 


Ax  =  {A,  H). 


In  fact,  Equation  (4.14)  is  recovered  if  A  is  replaced  by  u  and  v  in  the  above  equation: 


ux  =  {u,H}  =  —  = 

=  iu*  ve~1^1  +  iuyy 

Vi  =  {v-h}  =  W2  =  -'TOV^-n;, 

=  iu2e+l^x/2  +  ivyy.  • 


(4.23) 


In  addition. 


n2,  =  {nQ,//}  = 


yield  the  complex  conjugate  equations. 

Equation  (4.14)  may  be  recast  in  the  form  of  an  autonomous  system.  For  such  a 
system,  the  Hamiltonian  is  a  conserved  quatity.  Let  i;  =  ve~'^®x,  and  substitute  in 


Equation  (4.14),  resulting  in 


ux  -  ih\uyy  +  iu*v  =  0 


vx-iK2vyy  +  iAQv  +  iu2/2  =  0. 


(4.25) 


The  Hamiltonian  corresponding  to  Equation  (4.25)  is 


H  =  ft((ii’)2r]  +  A',K|2  +  K2\vy\2  +  ±Q\v\2 /2. 


(4.26) 
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where  9  reads  as  “the  imaginary  part  of”. 

4.1.2  An  Exactly  Solvable  Case 


(4.27) 


(4.28) 


(4.29) 


When  the  bottom  is  flat,  and  the  boundary  conditions  are  constant,  the  surface  system 
becomes 

alx  +  iKse~'^xa\a2  =  0 
a2l  +  iK%e+'^xa\  =  0 

(4.30) 

ai(z  =  0,y)  =  Ai 
a2(z  =  0,y)  =  A2 , 

where,  A3  are  constants.  The  above  system  is  very  familiar  to  the  nonlinear  optics 
community-c.f.  Arsmtronget  al.  [63].  Replacing  a,  =  /4,{x)exp0,(x)  in  Equation  (4.30) 
changes  the  system  to 

Ah  ~  A s-4i  A2  sin  ft  =  0 

(4.31) 

A2z  4-  KeA\  sin  ft  =  0 


HO 


for  the  real  part,  and 


Oh  +  A  $A2  cos  11 
•42#2x  +  KqA\  cos  n 
for  the  imaginary  part,  where  H  =  26\  -  9\  +  6x. 
imaginary  part, 


=  0 
=  0 


(4.32) 


Combining  the  equations  from  the 


fir  +  6  4*  [K^A\l A?  —  2A5A2]  cos  11  =  0. 


(4.33) 


Thus,  Equation  (4.30)  is  equivalent  to 


Ah  -  K5A\A2  sin  11  =  0 
.42r  +  A'6/1j  sin  H  =  0 
Hr  +  6  +  [h$A\/ A2  —  2A5j42]  cos  H  =  0 

4i(0)  =  A 
Ai(0)  =  Ai 

H(0)  =  Ho.  (4.34) 

To  continue,  let 

X  =  A2  sin  H  (4.35) 

T  =  .42cosH  (4.36) 

z  =  .4 (4.37) 

Multiplying  the  second  equation  of  Equation  (4.34)  by  sin  11,  using  the  third  expression 
of  Equation  (4.34)  it  can  be  discerned  that 


*r  =  -KeZ  -6Y  +  2KbY\ 


(4.38) 
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since  Xx  —  A21  sin  +  cos  fE  Multiplying  the  third  expression  of  (4.34)  by  .42  sin  fl. 

Yx  =  SX  -  2K5XY,  (4.39) 

after  making  use  of  the  second  equation  in  Equation  (4.34).  The  Z  variable  may  be 
eliminated  from  Equation  (4.38)  by  noting  that  from  the  real  part  of  the  original  system 
that  conservation  of  energy  is  given  by 

Ks{X2  +  Y2)  +  KeZ  =  Eq.  (4.40) 

Eliminating  Z  using  Equation  (4.40),  dividing  Equation  (4.38)  by  Equation  (4.39),  leads 
to 

dX  _  K5(X2  +  3T2)  -  SY  -  E0 

dY  (6-2KsY)X  ’  (  ' 

which  may  be  used  to  investigate  the  structure  of  the  phase  plane  of  Ai-  The  dynamics  of 
A\  follow  immediately  from  the  conservation  of  the  energy  constraint,  Equation  (4.40). 
Three  cases,  depending  on  the  size  of  the  detuning  parameter  6 ,  are  investigated.  A  plot 
of  the  detuning  parameter  as  a  function  of  frequency  and  0  is  shown  in  Figure  4.3,  for 
the  dispersion  relation  given  by  Equation  (2.46).  When  6/2y/Eo  ~  0,  the  phase  plane 
is  shown  in  Figure  4.4.  Note  that  dX/dY  =  0  and  X  —  0  gives  the  two  centers,  at 
(.V,  Y)  =  (0,  ±y/Eo/\/TKl).  Setting  Y  =  0,  dX/dY  =  0  gives  the  radius  of  the  bounding 
circle,  at  y/Eo/y/Ks,  beyond  which  the  orbits  diverge.  Additionally,  there  are  two  saddle 
points  at  (AT,  Y)  =  (±v/£o/\/A7)0).  Motion  along  the  limiting  circle  takes  place  in  such 
a  way  that  A\  =  0,  and  A2  =  Eo/s/Ks ■  If  we  start,  for  example,  with  A\  ^  0  and 
A2  =  0,  motion  in  the  plane  takes  place  along  the  line  Y  =  0  up  to  the  limiting  curve, 
the  phase  fi  is  then  equal  to  jt/2.  From  the  imaginary  part  of  the  original  system,  it  may¬ 
be  deduced  that  the  variation  of  fi  in  this  limit  is  described  by 


-  2K5Eq/2cosQ  =  0. 


(4.42) 
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Figure  4.3:  Detuning  parameter  dependence  on  u>i  and  13. 

The  transition  from  the  state  sin  ft  =  1  to  sin  ft  =  -1  occurs  along  the  limiting  circle. 
The  distance  x  at  which  this  transition  occurs  is  infinite,  but  it  can  be  estimated  by 
solving  Equation  (4.42).  The  solution  is 

ft  =  tan_1[exp(-2A'5£,y2z)  tanfto],  (4.43) 

and  hence  an  estimate  of  the  spatial  length  at  which  the  energy  of  the  first  harmonic 
makes  an  almost  complete  transition  to  the  second  harmonic  is 

L  as  1/2A'5£'q/2,  (4.44) 

which  shall  subsequently  be  seen  as  related  to  the  “interaction  length”.  The  variation  of 
the  amplitude  of  Ai  along  Y  =  0  may  be  discerned  from 


Mx  =  A'5^2  ~  £0. 


(4.45) 
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Figure  4.4:  Phase  plane  for  A2  for  S  ~  0.  X  axis  is  vertical. 

which  is  obtained  by  eliminating  Z  from  Equation  (4.38),  and  making  use  of  the  energy 
relationship. 

The  solution  of  Equation  (4.45)  is 

M  =  ( Eq/K5)1/ 2  tanh[A'51/2£o/2(x  -  x0)],  (4.46) 


with  A2  =  (Eo/h's)1/2  tanh[( K^EoY^xq).  At  the  begining  of  the  growth  process.  Aj  > 
A2  so  that  sinfi  =  1  and  the  growth  of  the  second  mode  is  independent  of  A2.  With 
the  solution  of  A2  in  hand,  using  Equations  (4.39)  and  the  first  expression  of  Equation 
(4.34),  it  may  be  shown  that 


Ai(x) 


_ Ai _ 

\J\  -  tanh^ATj^Eo^Xo] 


sech\ / KsEq{x  -  xq). 


(4.47) 


From  this  solution  it  is  concluded  that  irreversible  energy  conversion  takes  place  for 
<5  =  0.  This  solution  is  not  stable,  however,  since  the  stationary  states  are  reached  by 
motion  along  the  limiting  curve  on  the  phase  plane.  The  smallest  of  b  invariably  results 
in  motion  along  homoclinic  orbits  with  consequent  of  beat  in  the  amplitude  of  A\  and 
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1  /  2 

Figure  4.5:  Phase  plane  for  Ai  for  6/2 E0  <  1.  X  axis  is  vertical. 

only  one  is  possible,  and  the  energy  is  concentrated  mainly  in  the  lower  mode.  The  two 
modes  interact  weakly,  and  the  spatial  beats  get  smaller  and  shallower  as  the  detuning 
paramenter  is  increased.  In  fact,  when  6/2E0  >  1,  /li(x)  ~  A\,  and  >42(0)  =  the 
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first  two  terms  in  Equation  (4.33)  are  dominant,  so  that  the  phase  is 

Q  =  6x  +  7r/2.  (4.50) 

Substituting  the  above  expression  into  the  second  expression  of  Equation  (4.34)  we  obtain 

.42  =  A2  +  —~A\  sin  bx.  (4.51) 

V 

The  phase  portrait  for  this  case  is  shown  in  Figure  4.6. 


Figure  4.6:  Phase  plane  for  A2  for  6/2E »  1.  X  axis  is  vertical. 

Note  that  of  the  three  cases  considered  here,  the  only  ones  physically  relevant  io  the 
sandbar  generation  problem  are  the  first  two  cases.  The  large  detuning  parameter  case 
violates  assumptions  on  the  size  of  the  wavenumber/frequency  in  the  model.  However, 
much  is  to  be  learned  about  the  surface  system  from  looking  at  the  high  frequency  case 
in  detail. 

When  uj i  is  large,  or  equivalently,  when  b  is  large,  the  amount  of  energy  from  the  first 
mode  transferred  to  the  second  one  may  be  quite  small.  As  was  just  mentioned,  in  such 
case  the  first  mode  has  nearly  constant  amplitude.  Assume  that  the  boundary  conditions 


are  constant,  i.e.  0,(0)  =  Ay  Thus,  a\(x)  ~  <4t,  and  the  second  mode  expression  of 
Equation  (4.2)  may  be  integrated,  yielding 

<*2(2)  =  4^ie<5x-  (4.52) 

c 

Substituting  Equation  (4.52)  into  the  first  mode  equation 

aix(z)  ~  ie-^,  (4.53) 

which  can  readily  be  integrated  to  yield 


a  1(1)  =  A \eai ,  (4.54) 

where  a  =  K^K^/ b\A\\.  Thus  ai(x)  is  approximately  sinusoidal  with  a  wavelength 
proportional  to  |*4j|. 

For  a  nonzero  bottom  an  exact  solution  is  not  possible.  Consider,  however,  the  case 


AXx  -  A2sinfi  +  A\  =  0 


(4.55) 


■&2x  +  sin  fi  4-  X2-^2  —  0 

where  x:  represent  constants.  We  still  cannot  solve  this  system  analytically,  unless 
Xi  =  X'2  =  X,  in  which  case,  conservation  of  energy  is  given  by 


K5Al  +  K6A]  =  E0e'2xi. 


(4.56) 


Introducing  new  variables 

X  =  XElJ2e~xx 

Y  =  YEll2t~xx  (4.57) 


and  the  reduced  distance 
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assuming  6  =  0,  we  obtain,  using  Equation  (4.56)  the  system’s  phase  plane  equation 

dX  1  -  A'2  -  3F2 

-3v  =  XY  ■  (4'39» 


which  has  the  same  structure  in  the  phase  plane  as  that  shown  in  Figure  4.4.  The 
important  distinction  is  that  f  is  related  nonlinearly  to  x.  Whence,  the  damping  of  the 
waves  is  characterized  by 

9  =  x/(2A'5£o/2).  (4-60) 


For  x  "C  1»  there  is  weak  damping  and  the  waves  travel  a  considerable  distance  before  the 
energy  is  fully  dissipated.  On  the  other  hand,  if  \  >  1,  only  a  small  arc  of  the  trajectory 
in  phase  plane  is  traversed.  The  wave  substantially  attenuates  in  a  short  distance.  The 
relevant  case,  at  least  approximately,  to  the  oceanic  problem  considered  in  this  study 
is  the  former  case,  in  which  the  size  of  the  bottom  makes  the  coefficient  analogous  to 
X  in  the  above  presentation  of  0(a)  in  size  relative  to  the  other  terms  in  Equation  (4.55). 


With  an  understanding  of  the  dynamics  of  the  system,  we  now  present  the  analytical 
solution  to  this  special  case,  Equation  (4.30).  Our  development  follows  closely  Armstrong 
et  al.  [63].  For  the  sake  of  tidiness,  let  us  scale  Equations(4.30),  using 


w(x) 

v(x) 


l°i(*)l 

v/  K5E0 

M*)l 

V  K§Eq 

KsV  XsEqx 


AQ 


6 

KsEq 


(4.61) 


In  these  new  variables,  conservation  of  energy  assumes  the  simple  form 


v2(x)  +  w2(x)  =  1 


(4.62) 
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and  Equation  (4.30)  is  expressed  as 


dw 

—  =  —wvsinQ 

dx 

—  =  w2sinQ 

dx 

dfl  d  i 

—  =  i\Q  +  cot  fi  —  (ln(u>2t;)). 

dx  dx 


(4.63) 


Q  =  2&i (x)  -  02{x)  -f  A Qx  here.  For  A Q  =  0,  i.e.  perfect  phase  match,  Equation  (4.63) 
provides  another  constant  of  integration 


T  =  w2  vcosQ  =  tn2(0)t;(0)  cosfi(O). 


(4.64) 


Thus,  making  use  of  this  constant  of  integration,  and  conservation  of  energy,  it  readily 
follows  from  Equation  (4.63)  that 


^  =  ±2[).2(i-u2)2-r2]1/2, 

with  the  sign  being  determined  by  the  sign  of  sinSl( 0).  Hence 


(4.65) 


-±r 

Jv* 


v2(i) 


d(v)2 


(4.66) 


2(0)  [^(1  _  V2)2  _  r2]V2 

which  is  the  Complete  Elliptic  Integral.  Since  v  is  real  and  less  than  or  equal  to  1.  v2  is 
constrained  to  move  between  the  two  lowest  roots  of  v2(l  -  v2)2  -  F2.  Call  these  va  <  Vf,. 
We  then  arrive  at  a  general  definition  for  the  interaction  length-c.f.  Equation  (4.44)- 
which  is  the  spatial  expanse  in  which  the  solution  goes  from  one  root  to  the  other. 


r  _  _ d{v)2 _ 

*~Jvl  fu2(l  -  172)2  _  r2!1/2' 


(4.67) 


vl  [u2(l-t?2)2 

If  the  boundary  conditions  are  such  that  T  =  0,  v2  =  0,v2  =  1,  and  L ^  —  oc,  the 
solutions  would  be 


rr=0  =  tanh(x  +  x0) 
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u>r_0  =  sech(x  +  i0).  (4.68) 

The  case  x0  =  0  corresponds  to  A 2  =  0.  On  the  other  hand,  if  A\  =  0.  then  x0  —  oc. 

If  both  Aj  /  0  and  #2(0)  -  20i(O)  =  ±tt/2,  depending  on  sinQ(O)  =  ±1,  the  second 
harmonic  or  the  fundamental  gets  amplified  first.  If  the  fundamental  gets  amplified  first 
( io  <  0),  the  second  harmonic  decreases  to  zero  and  then  increases  until  all  the  energy 
is  in  the  second  harmonic.  If  the  second  harmonic  is  amplified  first  (£0  >  0),  there  is 
complete  energy  conversion. 


To  write  down  the  explicit  solutions,  define  vc  >  Vb  >  Vl-  the  third  root  of 


u2(  1  -  u2)2  -  r2  =  0. 


(4.69) 


Let 


j.2  =  O’2  ~  Va) 

(vi  -  viy 

72  = 


K2  -  vl) 


be  the  argument  and  modulus,  then 

±1  /**> 


V(vI  -  vi)  4(o)  [(l  -  y2)(i  -  72t/2)]1/2 

and  the  amplitude  squared  solutions  are,  in  terms  of  Jacobi  Elliptic  functions  “sn”, 


dy 


(4.70) 


(4.71) 


v2(x)  =  v2  +  (^-u2)sn2[(u2-u2)1/2(x  +  x0);7] 

w2(x)  =  1  -  v2(x),  (4.72) 


with  io  being  determined  by  the  boundary  condition  y2(0)  and  the  value  of  7.  Note  that 
Tmax  =  4/27,  as  determined  from  Equation  (4.72)  and  Equation  (4.69). 

The  solutions  for  imperfect  phase  matching,  AQ  ^  0,  can  be  found,  using  variation 


of  parameters.  Equation  (4.64)  is  now 
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the  first  case,  vi>  and  L\  are  relatively  insensitive  to  A Q  and  substantial  power  transfer 
occurs,  the  interaction  length  is  very  large.  On  the  other  hand,  for  A Q  >  1,  there  is  less 
power  transfer  and  the  interaction  length  is  shorter.  Figure  4.8  shows  how  the  interaction 
length  varies  nonlinearly  with  m  (and  hence  with  AQ),  and  in  Figure  4.9  and  Figure  4.10 
illustrate  the  dependence  of  the  interaction  length  on  the  size  of  the  nonlinear  parameter 
a,  and  the  dispersion  parameter  3-  The  relevant  size  of  the  parameters  a  and  3  in  the 
oceanic  application  under  consideration  is  as  high  as  0.15  for  a,  and  0.005  <  3  <  0.15. 
Hence,  from  the  graphs  it  may  be  inferred  that  the  interaction  length  is  more  sensitive 
to  dispersion  than  to  nonlinearity  for  the  above-mentioned  ranges  of  a  and  3- 


Figure  4.7:  v2  dependence  on  the  detuning  parameter.  In  all  cases  w2(x  =  0)  =  1.  The 
interaction  length  and  the  power  transferred  to  v2  decreases  as  A Q  increases. 


As  a  way  to  assess  the  evolution  of  waves  with  periodicity  in  the  longshore  direction. 
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Figure  4.8:  Interaction  length  dependence  on  the  nonlinear  parameter  AQ. 
suppose 

ai(x,y)  =  t i(x,y)e’(f!lI-u;il+,,y) 

a2  (x,y)  =  v(x,y)e'{k2X-U2t+hy)  (4.81) 

Then  the  system  (4.2)  is  now 

ux  +  il\Kxu  +  iA:5«*ue~,(7y+^)  =  0  (4.82) 

vx  +  iljK2v  +  iK6u2e+'Uy+Sx'>  =  0,  (4.83) 

where  7  =  /2  —  2/j,  which  can  be  zero.  Consideration  shall  be  made  here  to  the  high 
frequency  case.  For  u\  large,  u(x,y)  ~  u°.  Hence  Equation  (4.83)  may  be  integrated. 
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a 


Figure  4.9:  Interaction  length  dependence  on  the  nonlinear  parameter  a. 


yielding 


K6u  2e*(7v+^) 
6  +  ljK2 


(  1.84) 


assuming  v(x  =  0,  y)  =  0  and  u(x  =  0,  y)  =  u°  constant.  Using  this  expression  in 
Equation  (4.82), 


u(a:,t/)  =  u°exp[-z/i  A'ii  +  i 


K5K6\u°\2x, 


6  +  1* 


(4.85) 


Hence  v  oscillates  with  lines  of  constant  phase  normal  to  the  tan-1(^)  direction,  where 
the  angle  is  taken  with  respect  to  the  shoreward  direction.  When  I2  -  2 1\  exactly,  the 
direction  of  constant  phase  orthogonals  is  the  shoreward  direction.  On  the  other  hand, 
u  oscillates  in  the  x  direction,  with  a  repetition  length 


Ly  = 


6  +  l]K2 

2tr 


{A'sA'i 


61  “ 


-(^  +  /22A-2)/?A'1) 


-1 


(4.86) 


Furthermore,  v  can  develop  a  singularity  when 


6  +  l\  A  2  =  0, 


(4.87) 
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P 


Figure  4.10:  Interaction  length  dependence  on  thp  dispersion  parameter  /?. 


that  is,  when  12  =  K 2  (note  that  h  <  0).  In  terms  of  the  y  component  of  the 

wavenumber,  the  singularity  occurs  when 


h  =  ±1 


-26k7 


a 


(4.88) 


An  /2  of  such  value  is  not  at  all  unreasonable  to  consider.  A  singularity  must  be  in¬ 
vestigated  much  further,  since  it  is  most  likely  a  result  of  the  method  of  analysis  used 
here.  However,  there  is  a  change  in  sign  in  v  on  either  side  of  the  location  at  which  the 
singularity  is  predicted. 

Yet  another  interesting  feature  is  the  situation  when 


«°|2  =  (5  +  /22A'2)/^1/A'5A'6 


(4.89) 


again,  a  reasonable  value.  In  such  an  event,  the  modulation  of  u  practically  disappears. 
Then 

u(i,y)ss  u°  =  ±\J(6  +  l\K2)l\Kil A'5A'6  (4.90) 
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and 


v(x,y)  =  - 


l\Kxe'[bl+7y) 


f  4.9i : 


which  is  a  simple  sine  wave.  Hence,  one  could  conceivably  use  modulations  in  the  y 
direction  to  nonlinearly  produce  linear  sine  wave  signals  of  the  second  harmonic  with 
amplitude  given  by  l2.  If  l2  =  2 l\  exactly,  the  wave  oscillates  in  the  shoreward  direction. 

Carrying  out  this  high  frequency  analysis  further,  we  can  consider  the  effect  of  the 
bottom  topography  under  special  circumstances:  the  case  when  /(x,y)  =  f(y)  leads  in 
a  straight-forward  manner  to 

A'6u2e'(^y+^x) 


V  —  — 


6  +  i22K2  +  h\f(yy 


(4.92) 


(4.93) 


again,  assuming  v(x  =  0,  y)  =  0  and  u(x  =  0,  y)  =  u°  constant,  and 

,  ,  0  r  ,2r/  ,  •  A'5A'6|a°|2X  , 

Thus,  the  effect  of  the  bottom  in  this  case  is  to  change  the  amplitude  of  v,  and  at  the 
same  time  modulate  the  oscillations  of  u.  Again,  the  possibility  of  a  singularity  and  a 
change  in  sign  in  v  exists. 

Finally,  the  same  method  may  be  employed  to  assess  the  effect  of  a  mildly  sloping 
bottom  on  the  high  frequency  solution.  Assume  f(x,y)  =  vx/2,  where  v  is  small.  The 
same  procedure  leads  to 

(u0)2AV<'>'v+^) 


v  ~ 


[l-< 


6  +  K2ll 

2K<\v 


(S  +  Kil2)2 
u  =s  u°exp{-t(  A'i/j 
|u°|2A'5A'6 


{(*  +  Kill)2!2 /2  +  i(6  +  Kil2)x  -  i}]e-A>*2 
|u0|2A'4A'5A'6, 


)x  -  i(  A'3  - 


(6  + Kill)2'"  "“J  (6  +  K2I2)2 

2|  u°|  K^KsKeu  (6  +  K2l2)2x2 


exp  - 


(s  +  K2i2y 


■( X  - 


W} 
) 


(4.94) 
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The  result  is  only  valid  for  K4vx2  <C  1.  That  is,  since  A'4  is  of  the  same  order  as  .  it  is 
valid  for  x  «  0(  In  order  to  discern  what  is  fundamentally  different  about  the 

sloping  case,  consider  the  situation  in  which  u°  has  no  y  dependence,  so  that  Equation 
(4.94)  has  the  form 


(u°)2A'6e1^ 


[1-t 


2K4I/  b2x2 


b2  <y'  2 
,2  ,, 


-iSx-  l)]e-A'<1/r2 


0  r  -I  0 1'2  r-  r-  ■,  r-  \u°\2A4K5Ke,  2i 
u  exp{2|  uu|  K5K6x/S  -  i(h3 - - )vxz} 


Ill  u°[hi^K.P‘'(x-b\y2) 

e  b 


(4.95) 


From  Equation  (4.95)  it  is  readily  apparent  that  v  oscillates  proportionally  to  e'^r.  its 
maximum  amplitude  b/K6  times  smaller  than  u2.  The  phase  will  drift  quadratically  with 
distance  and  proportionally  to  K4v.  The  amplitude  drops  linearly  at  a  rate  proportional 
to  the  size  of  K4u  and  K3V,  the  wave  decays  exponentially  at  a  rate  controlled  by  the  last 
exponential  in  the  above  expression.  To  properly  interpret  the  decay,  recall  that  |<5|  >  1 
and  b  is  strictly  negative  in  this  analysis.  The  second  term  in  the  exponential  implies 
that  decay/blow  up  would  be  a  possible  outcome  of  the  original  model.  However,  this 
is  an  artifice  of  the  present  analysis.  If  the  assumption  u(x)  ~  constant  is  violated,  the 
above  expressions  are  not  valid.  Thus,  for  our  interpretation  to  be  valid,  it  is  required 
that  2|u|2A'4A'5A'6t'x/d3  <C  1. 


A  very  important  question  that  arises  in  the  applicability  of  slightly  resonant  inter¬ 
acting  triad  expansion  techniques  to  oceanic  waves,  is  that  we  may  be  neglecting  very 
important  side-band  modulations.  These  can  be  producing  interesting  structure,  control- 
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ling  the  stability  of  the  primary  waves,  or  affecting  very  minimally  the  structure  of  the 
evolving  waves.  A  general  result  on  this  issue  is  forthcoming,  but  for  now  we  limit  our 
attention  to  the  high  frequency  case.  The  problem  of  bands,  rather  than  isolated  modes, 
and  its  effect  on  the  evolution  of  individual  waves  has  been  examined  by  Hasselmann  [64] 
in  the  context  of  deep  oceanic  waves.  We  would  like  to  find  the  effect  of  side  bands  on 
the  main  waves  for  the  shallow  water  case.  The  following  analysis  follows  closely  work 
done  by  Brekhovskikh  and  Goncharov  [65]  in  connection  with  this  issue. 

Firstly,  the  modal  expansion  is  replaced  by  the  more  familiar  expression  for  the  lowest 
order  velocity 

/°C 

auJ(x)e~,UJtcLj  (4.96) 

•OO 

where  au)(x)  =  a(x.u>),  and  a'u(x)  =  a-u)(x)  since  u0  is  real.  Assume  a u  =  p^j  exp(  ik^x), 
where  ku  =  k(uj )  is  found  via  the  dispersion  relation.  Again  reality  means  that  k*.  = 
and  plj  =  Substituting  Equation  (4.96)  into  the  original  equations  and  using  the 

compatibility  conditions,  an  expression  for  the  amplitude  of  the  incident  waves  is 

d  f00 

S~Pu  =  -iotu)  /  pqps  e\p(-i^qs0Jx)dq  (4.97) 

ox  J — oo 

where  s  =  u  —  q,  and  A qsU)  =  q  +  s  -  k w.  If  the  incoming  harmonic  wave  u(0,<)  = 
cti(0)e-lU;i*  +  c.c.,  i.e.  Pu;(0)  =  ai(0)£(u;  -  uq)  4-  a\(0)6(u  -f  uq),  the  spectrum  of 
u(x,t)  remains  discrete  at  any  time:  the  only  non-zero  components  are  u;n  =  nj 
kn  =  fc(wn),  n  =  ±1,  ±2,  •  •  •  and 


au(x)  =  ^  an(x)«5(u;  -  u;n). 


Then  Equation  (4.97)  yields 


d 


u  V'  -t'Zj 

—an  =  -iaujn2_^a  ate 


mlnx 


(4.98) 


(4.99) 
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/  —  Tl  7TI,  Ajp/^  —  km  T  k[  hf±i  Cin  —  CLn ,  and  0.Q  —  0. 

Taking  u\  as  the  principal  harmonic  and  6lj  as  the  width  of  the  spectral  band,  we 
extend  Equation  (4.30)  to  include  the  interactions  of  spectral  components  of  the  wave 
train  with  long- wavelength  waves.  Except  for  a  constant  multiplying  the  integral,  the 
spectral  amplitude  equation  is 

d  f 

Ox  ^  MXUJ  I ^  [pu?i  +£PuJ—lJi  —  T  ]d£, 

A±  =  h^j  kuj ±  kuj^+^u)  ~  -f  kij1  ± 


±kuJ'  ±  ~LJ')  =  ku;  ~C3  u  =  lcp*  -  C9  ^  =  A“ 


(4.100) 


where  cg  is  the  group  velocity  and  cph  is  the  phase  velocity.  Approximating. 


Pu>i-K  ~  Puii+uj+i  ~  Pu> i 


P-iJi+\jJ+£  ~  P— U>]+£  ~  P— lx>i  ■ 


(4.101) 


the  equation  for  the  amplitude  is 


S~Pu  ~  -ia^lpui  |2etAu;xA uj. 
ox 


(4.102) 


As  w'as  done  in  the  discrete  case,  assume  the  frequency  is  sufficiently  high  so  that  p^j^ 
constant.  Thus, 


/V  u  I  1 2  Am  «|pw,  |2e,AwXAu2 

pu(x)  *  -a  — I  pUl  \2e  w  Au;  =  -  —  - 7T 

Au;  [c „fc  (0)  “  cg  Vl) 


(4.103) 


The  following  equation  for  corresponds  to  such  an  interaction: 


f  -iAfi  jf  «  |  pUl  |2pw(At*;)2 

pu  =  -mu  PiPu-(e  d£  *  -^-—7 

J  l\ul  lCp/,(0)  —  Cg  (u^l)j 


(4.104) 


where  (u;  -  u\  |  <  A^t.  Let  ai  =  pu^u  stand  for  the  amplitude  of  the  principal  harm- 


monic.  Then,  taking  into  account  the  term  corresponding  to  the  interactions  witli  the 
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second  harmonic,  we  obtain,  instead  of  Equation  (4.53). 

—a{  =  ia2uj2\a1  |2{A-1  +  -  c”1  (^.’i  )]_1  }ai -  (4.105) 

Its  solution  is  aj  =  a^OJe1171,  corresponding  to  waves  with  “spatial’’  shift  of  a  - 
-a2oj2\ a\  |2{il-1  +  w-t[cj/l1(0)  -  cj!(u;i )]-1}.  Hence,  in  the  high  frequency  the  main 
difference  between  the  discrete  and  the  banded  spectrum  case  is  that  the  latter  has  an 
additional  term  in  the  nonlinear  shift  as  compared  with  the  case  of  Equation  (4.53). 

4.2  Remarks 

The  conditions  for  the  stability  of  the  surface  system  and  of  the  full  model,  as  of  this 
writing,  have  not  been  analyzed  in  detail.  However,  it  is  possible  to  infer  from  the  results 
of  this  chapter  that  the  stability  of  the  surface  system  does  not  depend  on  the  frequency 
of  the  water  waves  since  only  weak  resonance  is  possible,  which  in  turn  means  that  less 
energy  is  shifted  from  the  lower  modes  to  the  higher  ones  the  higher  the  frequency  of  the 
fundamental  mode. 

Evidence  from  numerical  calculations  and  of  our  preliminary  analytical  work  on  the 
subject  suggests  that  the  stability  is  controlled  by  the  size  of  a,  by  the  possibility  of  a 
singularity  in  the  denominators  of  the  coefficients  Ks  and  A^6  by  the  right  combination 
of  parameters-see  Figures  5.3  and  5.4-  or  by  the  choice  of  boundary  conditions  Ai¬ 
ks  shall  be  shown  in  Chapter  6,  there  are  a  number  of  ways  in  which  refraction  occurs 
in  the  modes.  The  model’s  assumptions  places  a  restriction  on  the  degree  of  y  dependence 
of  the  solutions,  and  care  must  be  exercised  so  as  to  not  violate  the  assumption,  especially 
when  the  domain  involved  is  large.  It  may  be  possible,  however,  even  when  weak  y 
dependence  assumptions  are  not  violated,  for  the  solutions  to  lose  their  stability  due  to 
severe  refraction.  At  a  later  stage  in  this  study  we  shall  pursue  this  issue  with  the  hopes 


of  arriving  at  an  estimate  of  when  and  how  this  form  of  instability  occurs. 


Chapter  5 

Numerical  Solution  of  the  Model 


This  chapter  is  devoted  to  the  details  of  the  numerical  solution  of  the  full  system  and 
to  a  performance  evaluation  of  the  various  computational  schemes.  As  mentioned  previ¬ 
ously,  the  two  disparate  time  scales  effectively  decouple  the  fluid  and  sediment  problems, 
enabling  us  to  solve  the  full  model  iteratively.  The  input  to  the  model  is  comprised  of  an 
initial  bottom  configuration  and  the  mode  amplitudes  at  the  line  x  =  0.  The  required 
dynamical  parameters  are  the  fundamental  frequency,  an  estimate  of  the  size  of  a  and 
/?,  and  the  dimensions  of  the  rectangular  patch,  0  >  x  >  M,  0  >  y  >  N,  of  ocean  on 
which  the  solution  is  to  be  computed. 

Finite  difference  techniques  were  adopted  in  this  study  for  a  number  of  reasons,  the 
primary  one  being  that  they  are  very  well  suited  for  the  numerical  solution  of  the  hyper¬ 
bolic  equation  of  the  type  represented  by  the  mass  transport  equation.  Other  reasons 
have  to  do  with  practicality:  there  are  3  equations  to  implement  (the  mass  transport 
equation  plus  the  two-component  surface  system),  which  are  most  conveniently  solved  on 
the  same  grid.  In  addition,  synthetic  boundary  conditions  on  the  lateral  boundaries  were 
found  to  be  easily  handled  using  finite  difference  techniques.  In  this  study  a  uniform  grid 
was  found  to  be  adequate  for  our  purposes  and  hence,  used  exclusively. 
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5.1  Implementation  of  the  Mass  Transport  Equation 

The  mass  transport  equation  shall  be  implemented  numerically  using  a  Two-Step  Lax- 
Wendroff  scheme,  which  is  second  order  accurate  in  time  and  space.  Since  this  technique 
is  very  well  established,  we  shall  not  report  on  such  standard  issues  as  consistency, 
convergence,  and  uniqueness.  The  reader  is  directed  to  Smith  [66]  for  those  details. 

We  shall  define  the  following  difference  operators: 

A,  =  u(qJ+ 1)  -  u(qj)  forward  difference 

V,  =  u(qj)  -  n(m_i)  backward  difference 

(5.1) 

Sq  =  u(q]+ !/2)  -  u(qj_ j/2)  central  difference 
Aq  =  u(qJ+ 1)  +  u(qj)  forward  average 


in  the  independent  variable  q,  say.  The  physical  space  is  given  by  1Z2  x  T  =  [0  <  x  < 
A/,0  <  y  <  JV]x{T  >  0}.  Define  Tl\xT^  =  (xr,y3)xTn  =  (rAx.sAy)xnAT  €  TZ2xT. 
Furthermore,  there  are  integers  m  and  n,  such  that  M  =  mAx,  N  =  nAj/. 

Equation  (3.50)  is  approximated  by  the  following  computational  module: 


Ktl'/l  -  4^*  +  +  2AyAyl/ 

hl+1  =  hi  +  §7' 6*ti*  +  (5-2) 

on  72^  x  T^.  The  module  is  illustrated  in  Figure  5.1  for  one  space  dimension. 

For  the  sake  of  clarity,  the  stability  criteria  shall  be  established  in  the  shoreward 


direction  only.  Since  nx  =  Hhhx  substituting  in  Equation  (5.2)  yields 


y+l  =  h n  -  &[i(A,  +  Vx)fc  -  ^X(AX  -  Vx)ft] 


where  £x  =  -^Ar/2Ax.  A  local  stability  criteria  may  be  established  by  replacing 
h  =  £nexp(!rAx)  in  Equation  (5.3),  from  which  it  follows  that  the  growth  factor  is  such 
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Figure  5.1:  Computational  module  for  the  mass  transport  equation. 


that 

I  Cl2  =  l-Ul-^)d-cos(rAx))2.  (5.4) 


and  formal  linearized  stability  shall  result  if  jCI  <  1,  which  restricts  £2  <  1. 

Usinu  he  same  argument  the  stability  criterion  in  the  span-wise  direction  can  be 
found,  so  that  the  stability  of  Equation  (5.2)  in  two  space  dimensions  requires  that 
£  =  -(nhhxAT /2Ax,VhhyAT/2Ay)  be  less  than  one  in  component  form.  Since 


Ph. 

Vh 


2  4  32k3.  „  .  „  ~ 

-E- 7Ld-/32^;2/2)^|a;|2 

j=i 


A  4 32kJhl 


-  E  -“M1  "  ^h2kj/2)JAljR}y  -  IjyRjl 

]=i 


where  I}  and  Rj  are  respectively  the  real  and  imaginary  parts  of  ay  It  is  possible  to  show 
that  the  maximum  value  attained  by  either  |/'/,|2.  or  |  i>h  |2  is  of  the  order  of  1633|«;|4. 
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Hence,  for  stability  the  grid  size  is  determined  by  the  constraint 

f 


(5.6) 


a  result  which  sits  well  with  the  need  to  be  economical  about  computer  resources  and 
that,  as  will  be  referred  to  subsequently,  will  not  conflict  with  the  stability  criteria  of 
the  overall  iterative  scheme  of  the  full  model.  Thus,  in  component  form,  for  £x  <  1,  and 
assuming  |  a_,  |  <  1  over  the  whole  domain, 

AT 


<  yd“3a, 


Ax 


(5.7) 


and  for  £y  <  1,  the  same  argument  leads  to 


AT  a  3 

—  <  /T3. 

Ay 


(5.8) 


Dissipation  is  known  to  occur  except  when  =  1.  The  effect,  however,  can  be  quite 
small-fourth  order  in  Ax-if  the  wavelengths  are  restricted  to  being  much  greater  than 
the  grid  size  [66], 


5.2  The  Surface  Equations 

5.2.1  Numerical  Solution  of  the  Two-dimensional  Surface  System 

The  numerical  solution  of  the  two-dimensional  surface  system 


Qlr 

a2x 

ai(x  =  0) 
a2(x  =  0) 


-iA'3/(x)ai  -  ih'5e  ,*za*1a2 
-zA'4/(x)a2  -  iKee+'^za\ 
Ax 
A2, 


(5.9) 


where  A\  and  A2  are  constants,  shall  be  used  later  in  the  evaluation  of  the  imple- 
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mentation  of  the  three-dimensional  surface  system.  A  standard  fourth-order  explicit 
Runge-Kutta  scheme 

0r  +  \  =  Or  +  “(Pi  +  2P2  +  2P3  +  P4  ) 

0 

Pi  =  AxF(<t>r,Xr) 

P2  =  AlF(<J)r  +  iPl.Ir  + 

P3  =  AarF(<z>r  +  ^P2,xr  +  ^Ax) 

P4  =  AxF(<pr  4-  P3,xr  +  Ax),  (5.10) 

where  F  is  the  right  hand  side  of  Equation  (5.9),  and  the  vector  4>r  =  [ai(xr),a2(ir)]i  w^s 
used  to  numerically  solve  this  system.  For  details  on  the  applicability  of  such  a  scheme 
to  the  solution  of  Equation  (5.9)  we  refer  the  reader  to  Boczar-Karakiewicz,  Bona,  and 
Cohen  [lj. 

5.2.2  The  Three-Dimensional  Model 

For  the  surface  equations,  we  rewrite  Equation  (4.2) 

aXz  -  iK\a\yy  +  iKzf(x,y)a\  =  -iK$e~'^xa\a2 

a2x  -  iK2a2yy  +  ih\f(x,y)a2  = 

ai(z  =  0,y)  =  Ai(y) 
a2(x  =  0,y)  =  A2{y) 

(5.1 1 ) 

aiy(x,y  =  0)  =  0 
oiy(x,y  =  0)  =  0 
a2y{x,y=M)  =  0 
a2y{x,y=M)  =  0 

to  separate  the  linear  and  nonlinear  parts.  The  first  two  boundary  conditions  are  in¬ 
herent  in  the  physics  of  the  problem.  The  remaining  boundary  conditions  are  artificial. 
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These  Neumann  boundary  conditions,  combined  with  a  computational  procedure  that 
will  be  explained  presently,  ensures  that  the  overall  structure  of  the  solutions  shall  remain 
negligibly  affected  by  the  choice  of  lateral  boundary  conditions.  We  call  this  technique 
the  "zero-flux  procedure”. 

In  order  to  explain  why  this  procedure  is  needed,  let  us  spell  out  what  sort  of  problem 
we  are  faced  with:  Since  we  need  to  compute  a  solution  over  a  finite  domain,  care 
must  be  exercised  in  imposing  boundary  conditions  on  the  lateral  sides  so  that  we  avoid 
the  introduction  of  structure  in  the  solution  that  is  strictly  mathematical  rather  than 
physical1  in  nature.  To  avoid  this  situation  we  use  appropriate  boundary  conditions  along 
the  lateral  sides,  and  in  addition,  place  restrictions  on  the  initial  bottom  configuration 
and  the  boundary  condition  at  x  =  0  so  that  we  can  compute  an  oceanic  event  on  a 
swath  of  what  amounts  to  be  an  effectively  unbounded  domain.  We  have  found  that 
this  zero  flux  procedure  is  superior  to  other  synthetic  boundary  conditions  in  minimizing 
unwanted  structure  in  the  solutions. 

The  Neumann  boundary  conditions  make  the  problem  well-posed,  however,  by  them¬ 
selves,  would  introduce  a  great  deal  of  structure.  Physically,  these  boundary  conditions 
correspond  to  placing  hard  barriers  on  the  lateral  sides  of  the  domain.  A  posteriori  we 
know  that  the  solution  to  the  model  is  two  dimensional  if  neither  the  bottom  nor  the 
boundary  condition  at  x  —  0  has  y  dependence.  In  such  a  case  the  zero  flux  condition 
on  the  lateral  sides  has  no  effect  on  the  solution  over  any  part  of  the  domain,  i.e.,  it  does 
not  add  y  dependence  to  the  solution.  We  perform  the  calculation  of  the  system  over 
a  computational  domain  which  we  divide  into  three  regions.  The  large  central  region, 
flanked  by  two  sufficiently  wide  lateral  strips,  is  one  in  which  y  variation  in  the  initial 

1 A  possible  way  to  compute  a  solution  of  the  problem  over  an  effectively  unbounded  domain  over  a 
finite  grid  is  to  impose  periodic  boundary  conditions.  However,  periodicity  imposes  unwanted  symmetries 
on  the  structure  of  the  computed  solutions. 
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bottom  or  in  the  boundary  condition  at  x  =  0  is  possible.  In  the  lateral  strips  no  y 
dependence  in  the  above  mentioned  quantities  is  permitted.  The  solution  in  these  lateral 
strips  is  discarded.  The  initial  bottom  and  the  boundary  condition  at  x  =  0  are  con¬ 
nected  smoothly  in  all  three  regions  so  that  a  minimal  amount  of  structure  is  introduced 
in  the  solutions.  The  size  of  the  lateral  strips  is  determined  by  what  amounts  to  an 
educated  guess. 

An  efficient,  simple,  and  sufficiently  accurate  method  to  implement  the  above  non¬ 
stiff,  “locally”  nonlinear  system  numerically  is  now  the  focus  of  our  attention.  Several 
issues  have  motivated  the  particular  choice  of  scheme:  (1)  The  accuracy  requirements 
are  not  very  sophisticated,  since  we  wish  to  explore  a  phenomenological  question  rather 
than  achieve  engineering  accuracy;  (2)  a  uniform  grid  is  preferred  over  a  variable  one. 
so  that  both  the  surface  and  mass  transport  equations  may  be  easily  computed  on  the 
same  grid;  (3)  the  computational  domain  is  fairly  large  for  the  sort  of  problem  presented 
in  this  study.  The  method  presented  here  has,  among  its  best  features,  low  storage 
requirements  and  high  speed  as  measured  by  its  low  operation  count;  it  is  easy  to  imple¬ 
ment  on  conventional  hardware  using  recursive  data  structures,  and  to  a  certain  extent, 
parallelizable  on  vector  machines. 

Define  the  following  vectors,2  with  all  K's  real: 

k  =  i[h\,K2}T 

kf  =  if(x,y)[K3,  K4]t 

<t>  =  [a1(z,y),a2(x.y)]T  6  C2  (5.12) 

with  (x,y)  €  ft2^,  so  that  the  system.  Equation  (5.11),  is  now  recast  on  the  discrete  grid 


2The  superscript  T  means  “transpose”  in  what  follows. 
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Tl2^  as 


\dx  -  kdyy  +  kf]<z>  =  b (x,y,<p), 


(5.13) 


with  the  linear  part  on  the  left  hand  side  and  the  nonlinear  terms  on  the  right  of  the 
equals  sign,  plus  boundary  conditions, 


oy  =  0  on  y=0,  y=N 
<t>  =  d>0  on  x=0, 


(5.14) 


The  term  b(ar,d>)  represents  the  nonlinear  terms.  Succinctly,  the  above  equation  may  be 
written  as 

£d>  =  b  (5.15) 


where  C  is  the  linear  operator.  Let  L  be  a  suitable  discretization  of  such  linear  operator. 
Suppose  the  value  of  the  vector  <t>  at  level  r  for  all  s  is  known.  Making  use  of  fixed 
point  methods  the  value  of  the  vector  at  level  r  +  1  may  be  found.  Computationally,  the 
calculation  is  performed  in  two  steps:  let  /  be  the  index  of  the  iteration  and  let  4>  be  an 
intermediate  result.  Then  the  following  computational  scheme  is  proposed: 


L<p  =  b  (x,y,<t>1) 
L<t>l+X  =  b(x,y,<t>). 

Formally,  Equation  (5.16)  is  equivalent  to 


(5.16) 


L4>l+l  =  b(x,y,<t>1). 


(5.17) 


To  start  the  iteration,  the  value  of  the  field  variables  at  the  rth  level  in  x  is  used,  i.e., 
<t>°  =  <t>r ■  The  condition  for  convergence  of  Equation  (5.16)  is  ft  nd  by  appealing  to  the 
fixed  point  theorem. 
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For  the  purpose  of  determining  the  convergence  criteria,  define  C.  a  region  in  C 4.  the 
generalization  of  the  four  dimensional  real  space  to  complex  variables.  Let  $  and  h  £  C 
be  two  vectors  in  that  space.  Hence,  the  derivative  of  A  with  respect  to  <i>  is  defined  as 

a  dA, 

A*  =  J($)=  —  (5.18) 

If  the  second  derivative  is  continuous  for  all  $  £  C,  then  it  satisfies 

||  A**($,h.h)||<  A  ||  h  ||2  (5.19) 

for  all  $. 

Furthermore,  let  ||  •  ||p,  with  p  =  1,2,  oo.  represent  the  induced  norms 
L\  -  max1<J<n{j:f_1  Mol}  =  {£?=i  £"=  i 

(5.20) 

Loo  =  maxj<,<„{^"=1  Mol}. 

Finally,  define  a  super-system  on  1Z ^  as 

[dx  -  Kdyy  +  Kf]$  =  B(x,y,*)  (5.21) 

plus  boundary  conditions, 

=  0  on  v=0.  y=N 

(5.22) 

4>  =  *0  on  x=0 

composed  of  (5.11)  and  its  complex  conjugate,  with 

K  =  [k,k*jr  £  C 
Kf  (x,y)  =  [kf.kf*]r‘  ’ 

$  =  [a1(z,i/),a2(i,y),<(x,y),a5(x, y)}T  £  C.  (5.23) 


Let  L  be  the  resulting  discrete  operator  of  the  super-system,  composed  of  L  and  its 
complex  conjugate.  Choosing  L  non-singular  (hence  L  is  non-singular),  and  multiplying 
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both  sides  of  (5.21)  by  L  1 . 

$i+1  =  A  (x.y,$;). 

Define  the  iteration  discrepancy  as 

||  ^'+1  ||p=||  *'+1  -  *'  ||p  . 

Appealing  to  the  Fixed  Point  Theorem,  it  can  be  surmised  that 


(5.24) 


(5.25) 


P*‘+1||P  =  ||A(*1)-A(*,-1)||„ 

~  ||J(^,-1)^,||P<|| 

<  ||j(^-1)l|P||J(^-2)||P||^'||P<--- 

<  n  ||J(^)||P||^°||P,  (5.26) 

l=o 

provided 

0  <||  J(*')||P<  1.  (5.27) 

Equation  (5.27)  is  in  fact  the  convergence  criteria  for  the  iteration  process. 

To  establish  an  estimate  of  the  rate  of  convergence,  let  r  >  0  be  given  such  that  the 
set  of  vectors  S  =  {*:||*  -  s  ||p<  r}  contains  a  fixed  point  s  of  A(s),  i.e., 

s  =  lim  =  lim  A^1)  -  A(s).  (5.28) 

1  —  00  1  —  00 

Further,  let  «S  C  C,  J(s)  continuous  on  S  and  ||  J(s)  ||p<  1.  Then  there  exists  an 
s  >  0  such  that  the  fixed  point  iteration  is  convergent  whenever  ||  $°  -  s  ||p<  s.  Define 
||  e,+1  ||p,  the  measure  of  difference  between  the  (/  +  l)t/l  iterate  and  the  root.  Hence 

II  e'+1  ||p=||  *'+1  -  s  ||*||  J( s )e(  +  A"(C;e'.e')  ||p<||  J(s)e'  ||  +R  ||  e'  ||2  .  (5.29) 


Ill 


Quadratic  convergence  is  possible  if  J(s)  =  0, 


For  the  problem  in  question,  however,  the  best  rate  of  convergence  will  be  linear  since 


J(s)  ±  0: 


e'+1  |j» 


&Ti7r-s,JW 


A  measure  of  resources  required  in  the  computation  is  the  size  of  the  resulting  ma¬ 
trix  problem.  The  slightly  better  flexibility  in  the  choice  of  discretization  for  the  linear 
operator  L  is  the  key  advantage  of  this  method  over  others.  The  most  economical  dis¬ 
cretizations  are  those  that  lead  to  a  tri-diagonal  or  quinta-diagonal  matrix.  In  our  study 

Z/  =  (2AlAr_  2A?Vr)^-(A^^  +  kf*+l)<:>*+1’  (5-32) 

which  leads  to  an  n  x  n  tri-diagonal  matrix.  Its  computational  module  is  illustrated  in 
Figure  5.2.  and  is  commonly  known  as  the  Douglas  scheme.  L  has  eigenvalues 


r  1 


Figure  5.2:  Computational  module  for  the  linear  operator  of  the  surface  system. 
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A a  =  -(3  +  2 p  +  2Axkf )  +  2pcos[^j-]  s  =  1. •  ■  -n  , 


(5.33) 


and  the  eigenfunctions  are  given  by 


{sin 


SIT  .  2  Sir  .  STT  J 


- ,  sin - , 

n  +  1  n  +  1 


•  •  -  ,sin 


n  +  1 


}■*  s  =  1,  •  •  • .  n. 


(5.34) 


Furthermore,  the  operator  L  is  diagonally  dominant,  since 

n+1 

I iij  I  ^  I  f'nl  i  —  1 , •  •  •  n 
i*i 


(5.35) 


the  Lij' s  being  the  entries  of  the  matrix  L,  and  non-singular  since 


I  La  |  >  |  La+i  |  >  0  i  —  1,  •  •  1 

I  La  |  >  |  L  |  +  |  L{t- 1 1  I„+i  Z„_i  ^  0  *  =  2,  •  •  -  n  —  1 

|Xr:|  >  i  =  2,---n.  (5.36) 

If  0  —  fre‘*,  where  0  =  aAys,  and  p  =  2^^-k,  upon  substituting  these  quantities  in  L 
the  magnification  factor  is 

5  =  2,(1  -  cos«)  +  2A*kr  +  3{2  *  V/l  -  W  -  -  2Alkf)’  <5  37> 

from  which  it  is  clear  that  |£|  <  1.  Thus  the  linear  operator  is  unconditionally  stable. 

An  estimate  of  the  accuracy  of  the  discretization  of  the  linear  operator,  and  a  check 
on  its  consistency  with  the  continuous  operator  on  the  grid,  is  given  by 

Az2  Ay2 

(L-C)<t>=  —  d>xr +  k-^-dW +  •••,  (5.38) 

where  C  is  the  continuous  linear  operator.  Equation  (5.38)  implies  that  the  scheme  is 
0(Az2  +  A y2)  accurate. 

This  order  of  accuracy  is  an  upper-bound  on  the  accuracy  of  the  overall  scheme,  hence 
attempting  to  reduce  the  error  ||  e1  |[  much  below  this  is  pointless.  Since  the  real  root  is 
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not  known  a  priori ,  the  iteration  procedure  is  carried  out  until  we  are  safely  below  the 
above-quoted  error,  but  not  much  beyond  that.  This,  in  effect,  is  the  criteria  used  in  the 
code  for  stopping  the  root-finding  iteration  procedure. 

Consistency  of  the  discretization  is  readily  established  by  comparing  the  continuous 
problem  with  its  discretization  in  the  limit  as  the  grid  size  gets  smaller.  It  may  be  shown 
that  the  discretization  approaches  the  continuous  operator  on  the  grid  uniformly. 

Since  the  Douglas  scheme  is  inapplicable  at  r  =  0,  a  standard  Backwards  Euler 
scheme, 

—  Ay?  ^  ^■r(kfd’)r-l-l ,  (5.39) 

is  used  to  discretize  L  for  the  first  step  in  x,  which  can  be  shown  to  be  unconditionally 
stable  as  well. 

Having  made  a  choice  on  the  particular  form  of  the  operator  L,  the  condition  that 
||  J(x)  ||p<  1  for  the  surface  system,  must  be  determined  explicitly,  so  that  convergence 
is  established  for  the  sand  ridge  problem.  To  estimate  the  size  of  J(x)  wc  make  use  of 
the  super-system,  Equation  (5.24), 


$,+1  =  L~lB(*')  (5.40) 

(5.41) 

$/-1  =  etc.  (5.42) 

Thus, 

<5$i+1  =3  J(*')6$'  (5.43) 

a  (5.44) 

etc.  (5.45) 
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with 

J  =  (5.46) 

where 

0  -iKhe-tf,XT+'a\l  -iK5e-i6xr+'al2  0 

-i2Kee-l6x'+'a[  0  0  0 

iK$e+iSxr*la!f  0  0  iK5e+i6x'+'a[ 

0  0  i2K6e+i6xr+'al{  0 

(5.47) 

for  the  Ith  iterate.  In  Equation  (5.47),  it  is  understood  that  the  a's  are  only  defined  on 
the  grid. 

From  Equation  (5.46),  ||  J  ||p<  1  if  the  size  of  L  is  greater  than  the  size  of  B'.  In  the 
L 2  norm,  the  convergence  condition  is 

II  •>  IIj=II  L-'B'  ||2<||  IT1  ||2||  B'  |b<  1,  (5.48) 

but 

II  L"1  ||2<||  I  lb  /  II  L  ||a=  2/  ||  L  ||2  .  (5.49) 

Replacing  (5.49)  into  (5.48)  yields 

II  J  ||2<  2  H  B'  ||2  /  ||  L  ||2<  1.  (5.50) 

Working  out  the  above  inequality  gives  the  condition  for  convergence  in  our  particular 
case: 

2^(Kl  +  4KZ)\a1\2  +  Kj\a2\2/  ||  L  ||2<  1.  (5.51) 

Since  LL*  =  L+L,  where  is  the  Hermitian  matrix  of  L,  then 


(5.52) 
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or  using  Eq.(5.33), 

v/(A'52  +  4A62)|a1|2  +  A-52ja2|2/(3  +  4 p  +  2Axkf)  <  1.  (5.53) 

where  the  L ^  in  y  is  used  to  estimate  the  size  of  the  vectors,  i.e.,  a,  =  maxi<,<„  rtf, 
i=l,2.  Hence,  Equation  (5.53)  gives  constraints  on  p ,  Ax,  and  a,,  to  be  satisfied  in  order 
to  guarrantee  convergence  in  the  solution  of  our  problem.  Another  constraint  we  impose 
in  the  numerical  implementation  is  to  restrict  Ax  to  be  less  than  2x/6,  so  as  to  minimize 
the  phase  error.  In  Figure  5.3  and  5.4  the  parametric  plots  of  A'5  and  K&  are  shown,  and 
are  included  to  complete  the  picture  of  the  relevant  size  of  all  the  quantities  involved  in 
Equation  (5.53). 


Figure  5.3:  Plot  of  A’s/a  versus  the  fundamental  frequency  uq,  and  Q. 

It  must  be  noted  that  owing  to  the  nature  of  the  nonlinearity  of  our  particular  prob¬ 
lem,  we  had  to  rely  on  the  super-system  to  arrive  at  a  convergence  criteria,  but  we  do 
not  actually  use  it  for  the  computation  of  the  field  variables.  Note  also  that  although 
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Figure  5.4:  Plot  of  K^/a  versus  the  fundamental  frequency  wj,  and  /3. 

inverse  matrices  appear  establishing  the  estimates  of  convergence,  they  are  not  actually 
required  in  the  computation  of  the  field  variables  themselves. 

5.3  Performance  Evaluation  of  the  Numerical  Schemes 

5.3.1  Evaluation  of  the  Mass  Transport  Equation  Scheme 

For  the  Lax-Wendroff  Scheme,  we  ran  a  few  test  runs  in  order  to  confirm  qualitatively 
the  scheme’s  stability,  consistency  and  accuracy,  checking  for  agreement  with  the  well 
established  theoretical  results.  Of  more  concern  to  us  was  the  issue  of  damping  and  of 
phase  drift.  In  order  to  quantify  the  scheme’s  dissipation  and  drift  we  used  a  model 
problem  for  which  an  exact  solution  is  known. 

The  model  problem  used  was 


hj  +  khhx 


0,  x  e  TV,T  >  0 


(5.54) 
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with  initial  condition 


1  i  <  0 

h(x,0)  =  <H-£i  0  <  x  <  l 


(5.55) 


y  1  +  si  X  >  l 

in  which  0  <  k  <  1,  and  £  <  1.  The  exact  solution  of  Equations  (5.54)  and  (5.55)  is 


x  <  kT 


h(x,T)=l  1  +  £fmr  kT  <x  <l  +  k(l  +  el)T 


(5.56) 


1  +el 


otherwise. 


We  tried  different  values  of  k  -  it  scales  the  time  step-,  but  we  report  our  results  for 
k  =  0.1,  and  for  such  a  case  convergence  is  possible  if  hAT /Ax  <  10  in  the  time  interval 
0  to  T.  Since  Equation  (5.54)  conserves  a  quantity  proportional  to  h p,  we  compared  the 
computed  value  h A  against  the  theoretical  value  h  as  a  function  of  a  =  kAT/Ax  and  as 
a  function  of  time  T  to  get  an  idea  of  the  scheme's  dissipation.  Specifically,  we  monitored 
the  constant  of  motion 

M/Ax 

c(T,a)=  £  h2A(T,rAx)rAx+±kT[h\(T,M)-h3A(T,0)},  (5.57) 

r=0  6 

where  M  is  a  very  large  value  in  xT.  For  an  estimation  of  the  phase  drift,  we  computed 


e2(a,  T)  =  £  |  hA(T ,  xr)  -  h(T ,  xT)  |2/  £  |  h(T,  xT)\2. 


(5.58) 


Figure  5.5  and  5.6  show  parametric  plots  of  c  and  e2  respectively. 

5.3.2  Performance  of  the  Runge-Kutta  Scheme 


The  accuracy  and  dissipation  of  the  explicit  fourth-order  Runge-Kutta  was  investigated. 
The  domain  was  128  units  in  extent,  or  roughly  10  interaction  lengths.  To  measure 
the  dissipation,  we  monitored  the  energy,  given  by  Equation  (4.62).  This  quantity  was 
conserved  by  all  trials  to  within  2%  for  all  reasonable  grid  sizes. 
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timeT 


Figure  5.5:  Dissipation  as  a  function  of  a  and  T  with  k  =  0.1  for  the  Lax-Wendroff 
Scheme.  From  top  to  bottom,  a  =  0.4,  0.2,  0.1,  0.05  respectively. 


To  estimate  the  accuracy  and  error  of  the  scheme  we  compared  the  outcome  of  the 
numerical  solution  to  the  exact  solution  (4.72)  using  the  following  measures: 


max{Erlx(rAJ)-  x'(*r)|} 

max{£r  |x'(rAx)|} 

Hr  I  x(rAx)  —  x'(xr)| 
£rlx'(rAx)| 

Erlx(rAx)-/(xr)|}2]1/2 

Er|x'(rAx)|2p/2 


(5.59) 


where  \  >s  the  calculated  value  of  a,,  and  the  exact  value  at  the  grid  location.  The 
exact  solution  was  computed  using  the  algorithm  given  in  [67],  pl89.  The  error  as 
a  function  of  grid  size  is  plotted  in  Figure  5.7,  from  which  one  can  conclude  that  the 
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timeT 


Figure  5.6:  Phase  drift  for  the  Lax-Wendroff  scheme  as  a  function  of  a  and  T  with 
Ar  =  0.1.  From  top  to  bottom,  a  =  0.4,  0.2,  0.1,  0.05  respectively. 

scheme  is  in  fact  fourth  order  accurate  and  consistent,  i.e.  the  error  drops  more  or  less  by 
a  factor  of  24  every  time  the  grid  size  is  halved.  For  the  accuracy  and  dissipation  trials 
=  0.5,  Ai  =  0,  in  Equation  (5.9),  a  flat  bottom  and  parameters  a  =  0.3,  /3  =  0.1, 
uq  =  0.5,  were  used.  Roughly,  10  interaction  lengths  was  the  extent  of  the  domain. 

5.3.3  Fixed  Point  Method  Performance  and  Evaluation 

Since  we  do  not  have  an  exact  solution  to  the  three-dimensional  system  we  sought  to 
discern  the  accuracy  of  the  Fixed  Point  Method  (FPM)  using  local  analysis.  Let  A  be 
the  size  in  x  or  y  of  each  grid  element-as  mentioned  previously,  the  grid  size  is  uniform 
in  the  domain.  A  comparison  of  the  computed  solution  at  a  particular  point,  using  A. 
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Figure  5.7:  Error  as  a  function  of  grid  size  for  the  Runge-Kutta  method.  l\:  . 

1 2:  —  —  —  —  ,  /qq  :  . 

and  a  solution  with  grid  size  A/2  yields 

I  *A  "  XA/2 1  s  *1  =  CO,[(A/2)"].  (5.60) 

Halving  the  grid  size  again 

IXA/2  -  XA/4I  =  =  CO[(A/*n  (5.61) 

Thus,  using  Equation  (5.60)  and  Equation  (5.61)  we  can  solve  for  p  to  get  an  estimate 
of  the  order  of  accuracy  of  the  scheme: 

log  k\  log  k2 

P~  log  2 


(5.62) 
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Using  the  same  parameters  and  boundary  conditions  as  those  used  in  connection  with 
the  Runge-Kutta  scheme,  and  a  domain  with  length  of  128  and  span  of  32.  we  found  that 
the  Fixed  Point  Method  _ '  ids  an  average  value  of  p  =  1.8,  with  a  standard  deviation  of 
0.5.  Values  of  both  fieM  quantities  were  used,  and  they  were  taken  from  various  points 
in  the  domain. 

We  did  not  perform  a  systematic  study  of  the  convergence  of  the  method.  We  ob¬ 
served.  however,  that  the  computed  values  tended  to  converge  as  the  grid  size  was  re¬ 
fined.  Since  comparisons  of  the  computed  solutions  with  an  exact  expression  for  the 
three-dimensional  case  were  not  possible,  we  compared  the  cross-sectional  values  of  an 
effectively  two-dimensional  solution  computed  using  FPM  along  the  whole  length  in  x 
and  midway  in  the  span-wise  direction  y,  with  a  solution  computed  using  the  Runge- 
Kutta  Method  with  a  very  fine  grid  spacing.  A  measure  of  the  error  is  given  by  the 


,  _  max{£r  k(rAx,muf)-  \7(x^) ]} 

°°(  '  V)  max{ErUVAx)|} 

MAX,  A,)  =  E.IvnAx.^-Vlx.ll 

lx  (rAx)| 

,  ,  ,  Er  I XirAx.mid)  -  \/(arr ) l}2]1^2 

=  E,l«va*)|’r 


MA  .Ay)  = 


(5.63) 


norms,  where  \  represents  the  solution  obtained  using  FPM  and  \'  the  solution  computed 
with  the  Runge-Kutta  scheme. 

For  the  case  Ax  =  Ay,  the  result  is  shown  in  Figure  5.8.  The  same  result  is  obtained 
when  we  fix  Ay  =  0.25,  and  we  vary  Ax.  On  the  other  hand,  when  Ax  =  0.25  is  fixed 
and  Ay  is  varied,  very  little  sensitivity  in  the  norms  is  obtained.  In  this  last  case,  the 
l\  4.4  x  10"  \  l\  4.4  x  10"3  .  and  %  4.4  x  10"3  for  all  sizes  of  the  grid  in  the  y 
direction  that  we  tried  3. 

The  rate  at  which  the  iteration  converges  in  FPM  as  a  function  of  the  grid  size  was 


’Note  that  there  is  no  y  dependence  in  the  solution  for  this  particular  trial. 
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Figure  5.8:  Error  as  a  function  of  grid  size,  with  Ax  =  Ay.  - -  -  lx: . , 

h- - ■ 

also  investigated.  With  a  =  0.3,  (3  =  0.08,  u>\  =  0.5,  and  boundary  conditions  A\  =  0.5 
and  Ai  —  0.1,  and  a  flat  bed,  we  monitored  the  iteration  discrepancy  at  a  particular 
value  of  x  in  a  fairly  large  domain.  As  expected,  we  found  that  the  number  of  iterations 
required  to  meet  a  certain  iteration  tolerance  decreases  as  the  grid  was  refined.  Figure 
5.9  shows  how  the  quantity 

n 

log10[max{^  ] 0l+l{x,sAy)  -  <t>l(x,sAy)\}}  (5.64) 

5  =  0 

drops  after  each  iteration  l  for  a  number  of  different  grid  sizes.  It  is  evident  from  the 
graphs  that  a  finite  and  small  number  of  iterations  are  required  to  reach  adequate  error 
tolerances  using  reasonably-sized  grids. 
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We  examined  the  iteration  convergence  at  the  first  step  in  x  as  well.  Recall  that  for 
the  first  step  a  Backwards  Euler  scheme  was  used  to  discretize  the  linear  operator  instead 
of  the  Douglas  scheme.  We  found  that  the  number  of  iterations  was  roughly  double  the 
number  required  elsewhere  in  the  domain,  where  the  Douglas  scheme  is  used. 
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Figure  5.9:  Iteration  discrepancy  as  a  function  of  grid  spacing.  The  number  of  iterations 
drops  as  A  =  4,  2,  1,  0.5  respectively. 

Before  we  examine  the  model's  speed  and  storage  requirements,  we  shall  present 
an  overview  of  the  implementation  of  the  surface  equations  using  Newton's  Iterative 
Method,  the  point  being  that  a  comparison  of  the  standard  approach  with  FPM  enables 
us  to  make  specific  claims  regarding  the  resource  economy  of  the  Fixed  Point  Method. 

In  the  most  straightforward  application  of  Newton's  Method  we  either  use  the  super¬ 
system.  or  separate  the  regular  system  into  real  and  imaginary  parts.  We  shall  opt  for 
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the  second  alternative.  Let 


G]  =  u  +  i  v 

a  2  =  w  +  iz  (5.65) 

and  use  some  suitable  discretization,  such  as  the  Douglas  or  Backwards  Euler  scheme, 
say.  Let 

F  =  {f\y  h,  h,  f\)  =  0  (5.66) 

represent  the  four  resulting  equations -here  we  have  placed  the  nonlinear  terms  of  Equa¬ 
tion  (5.11)  on  the  left  hand  side  of  the  equals  sign-for  the  values  of  the  field  variables 
at  level  r  +  1.  If  a  second  order  implicit  discretization  of  the  operator  dyy  is  used,  such 
as  would  be  the  case  if  the  Douglas  or  Backward  Euler  schemes  were  implemented,  the 
vector  F  has  the  following  dependence: 

F  =  F(xr+1),  (5.67) 

with 


X'  =  (U3'\US,  U3+l;V3-\V\  V3+l-,W3-\W3,  W3+l\  Z3~\ Z3,  Zs+l  )l  (5.70) 

is  the  ltk  iteration  to  the  approximation  X  of  the  exact  solution  x  which  is  being  sought, 
and 


*x'  =  (£s,^,75,n- 


(5.71) 
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Expand  F  about  X*  to  Linear  terms. 

F(xr+1)»  F(X')  +  J(X')-AX'.  (5.72) 

The  resulting  linear  system  for  the  unknown  vector  bXl  is  to  be  solved.  Let  X/+i  = 
X'  +  bXl  and  iterate  until 

maxi<s<n  ||  6Xl  ||<  0(f[Ax.  Ay))  (5.73) 

where  0(f(Ax.  Ay))  is  the  order  of  accuracy  of  the  discretization.  As  a  first  guess,  we 
set  X°  =  xr. 

The  Linear  system  resulting  from  an  implicit  scheme,  which  needs  to  be  solved  and 
recomputed  after  each  iteration,  has  the  following  structure: 

[A„ D.HX'  =  [/>,].  t=l,4  (5.74) 

where  for  each  i,  one  of  the  matrices  in  A  through  D  is  tri-diagonal,  and  the  other  three 
are  diagonal.  Hence,  the  full  system 

A/x  =  b  (5.75) 

is  4rc  x  4n,  and  while  sparse,  it  has  a  large  bandwidth. 

If.  on  the  other  hand,  an  explicit  scheme  was  used,  the  matrix  to  be  computed  and 
solved  for  each  iteration  would  have  been  4x4  in  size  and  full.  In  such  a  case  it  was  found 
that  the  accuracy  is  only  first  order  in  x.  and  the  grid  spacing  must  be  impracticably 
small. 

In  order  to  compare  the  Fixed  Point  and  the  Newton  schemes  insofar  as  economy  of 
resources,  we  need  to  present  details  of  the  solution  of  the  matrix  problem  in  Equation 
(5.75).  A  way  to  efficiently  solve  Equation  (5.75),  which  by  no  means  is  implied  to  be 
the  optimum  way,  is  to  use  a  pre-conditioning  matrix.  See  [68],  p527. 


126 


For  each  iteration,  we  need  to  solve  Equation  (5.75),  which  can  be  recast  as 


Mx  =  b. 

in  which 

M  =  C~XMC~ 1 
x  =  Cx 
b  =  C-1b. 


(5.76) 


where  C  is  a  pre-conditioning  matrix  such  that 

•  ||  C~x  M  ||<||  M  ||,  for  convergence, 

•  The  condition  number  k(C~xM)  <  k(M),  where  k(M)  =||  M  ||  /  ||  M~x  ||, 

•  C  is  easily  invertible, 

•  C  optimally  has  small  storage  requirements. 

The  iterative  matrix  solution  process  itself  is  thus 

x"+1  =  x"  +  u(C-lb  —  C_1Mx<).  (5.77) 

The  two  goals,  which  can  be  in  many  instances  incompatible,  are  high  speed,  measured 
in  number  of  computations,  and  low  storage  requirements.  In  what  follows,  we  shall 
compromise  on  storage  economy  for  the  sake  of  speed,  i.e.  suppose  the  computational 
domain,  which  is  always  fairly  large,  is  not  too  large.  A  good  choice  for  C,  since  M  is 
strongly  diagonally  dominant  in  our  problem,  is  to  use  the  symmetric  positive  definite 
tri-diagonal  part  of  M.  In  order  to  achieve  efficiency,  the  key  is  to  juditiously  carry  out 
the  multiplies  of  Equation  (5.77),  so  that  operations  are  performed  on  vectors  as  soon 
as  possible,  rather  than  matrices.  For  example,  after  computing  C-1,  which  incidentally 
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will  be  a  full  matrix  but  that  needs  to  be  computed  only  once  for  each  xr,  we  find 
Mxl  =  q,  a  vector;  next,  compute  C_1q,  another  vector,  and  so  on. 

The  operation  count,  for  the  Newton  Method  can  be  estimated  as  follows:  for  each 
value  of  xr,  we  need  to  solve  iterably  for  the  Ith  value  of  bX,  and  for  each  6Xl  we  need 
to  iterably  solve  Equation  (5.77).  First,  since  M  is  a  4n  x  An  matrix,  we  require  at 
most  0(n 2)  operations  to  find  C~1.  This  calculation  needs  to  be  done  only  once  for  the 
solution  of  Equation  (5.77)  since  we  use  the  same  C-1  until  the  iteration  of  Equation 
(5.77)  converges,  for  each  i/,  we  have  0(n2)  operations.  This  number  of  operations  is 
in  turn  performed  v  times  to  find  the  (v  +  l)th  vector.  The  estimate  for  the  number 
of  iterates  needed  to  solve  for  x"+1  depends  on  the  specific  form  of  M.  However,  we 
guess  that  the  number  of  iterations  required  is  perhaps  as  good  as  a  conjugate  gradient 
method,  which  is  typically  of  O(n).  In  addition,  we  need  to  iterate  /  times  to  get  the 
/  +  1  iterate  of  6Xl.  An  estimate  of  this  number  is  hard  to  estimate,  but  we  expect  that 
the  Newton  converges  quadratically,  whereas  the  Fixed  Point  Method,  as  we  showed 
previously,  converges  linearly.  Finally,  we  require  this  whole  process  be  performed  at  all 
values  of  x,  m  times.  Hence,  all  told,  we  have  mxuxlx  ( 0(n2)  +  0(n ))  operations,  and 
if  we  suppose  ^  =  O(n),  we  conclude  that  the  total  count  is  approximately  mxlx  0(n 3). 

An  estimate  for  the  operation  count  for  the  FPM  is  as  follows.  Equation  (5.16)  leads 
to  the  problem 

Lq  =  b  (5.78) 

for  the  unknown  d>,  where  L  is  a  2 n  x  2 n  tri-diagonal  matrix,  m  times  to  cover  all  values 
of  x  in  the  domain.  The  efficient  way  to  solve  Equation  (5.78)  is  to  decompose  the 
problem  in  two  steps:  let  L  —  WV,  where  W  is  a  lower  triangular  matrix  and. V  and 
upper  triangular  matrix.  Then,  solve 


ITg  =  b 


(5.79) 
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is  solved  for  g,  followed  by 

Uo  =  g.  (5.80) 

to  finally  obtain  0.  The  total  operation  count  for  the  solution  of  Equations(5.79)  and 
(5.80)  is  (5n  —  4)  multiplies  and  (3n  -  3)  adds.  AD  told.  0(l6n)  operations.  In  turn, 
this  process  is  performed  /  times  to  compute  the  {l  +  l)th  iterate,  and  finaUy  m  times 
to  cover  all  values  of  x.  The  total  is  m  x  l  x  0{n).  Thus,  the  operation  count  ratio  for 
these  two  methods  is  0(\/l  x  n2),  having  assumed  that  the  Newton  method  converges 
quadraticaDy  in  the  iteration  process.  Hence,  the  Fixed  Point  Method  is  considerably 
faster. 

We  can  also  compare  storage  requirements.  For  FPM,  we  need  to  store  the  old  and 
the  new  vector  at  each  x,  and  another  vector  for  the  iteration  process,  hence  we  store  6n 
values-note  that  for  our  problem  each  entry  in  the  vector  is  of  size  2n.  In  addition,  we 
need  to  store  all  the  entries  of  a  tri-diagonal  matrix  of  size  2n  x  2 n,  or  roughly  6n  values.4 
The  total  is  thus  12n,  or  O(n)  values.  For  the  Newton  Method,  we  have  12 n  values  which 
constitute  the  old  and  new  vector  elements  plus  a  storage  vector  for  the  iteration  plus 
roughly  '24n  for  the  matrix  M  entries,  with  an  additional  n2  for  the  pre-conditioning 
matrix.  The  total  is  0(n2).  The  ratio  of  storage  requirements  of  the  two  methods  is 
0(n). 

The  FPM,  as  we  hope  we  have  been  able  to  show,  has  many  attractive  features. 
Note  that  its  economy  of  resources  hinges  upon  the  simpUcity  of  the  matrix  that  the 
discretization  generated.  If  higher  order  accuracy  is  required,  the  matrix  will  probably 
be  more  compDcated  than  the  simple  tri-diagonal  matrix  that  was  used  in  this  study, 
requiring  greater  computational  resources.  A  somewhat  unavoidable  problem  with  the 

4 In  fact,  we  could  be  even  more  economical  and  use  multipliers  in  the  entries  of  L,  so  that  onl\  one 
half  of  the  tri-diagonal  matrix  entries  need  to  be  stored. 
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FPM  is  that  the  method  has  noticeable  dissipation.  As  we  shall  show,  however,  the 
dissipation  can  be  made  tolerable  at  the  expense  of  greater  computational  resources.  We 
do  not  know  with  certainty  the  cause  for  the  dissipation  and  further  study  is  required, 
perhaps  by  applying  this  scheme  on  a  nonlinear  equation  for  which  a  great  deal  more  is 
known  about  its  behavior  and  its  solution. 

To  illustrate  the  degree  of  dissipation  in  the  surface  system  FPM  implementation  we 
used  the  same  parameters  and  domain  that  was  used  in  connection  with  the  iteration 
issue,  and  we  fixed  the  iteration  tolerance  at  10-6.  Two  types  of  trials  were  carried  out. 
both  were  carried  out  using  a  flat  bottom.  In  the  three-dimensional  trial  we  assumed  the 
boundary  conditions  were  A\  =  0.5  +  0.0 ly,  and  Aj  —  0.1  +  O.Oly,  and  monitored  the 
conserved  quantity.  Equation  (4.26),  along  the  length  in  the  x  direction,  midway  in  the 
span-wise  direction.  The  derivatives  that  appear  in  Equation  (4.26)  were  second-order 
center-differenced.  In  the  two-dimensional  trial,  we  set  A\  -  0.5,  and  A2  =  0.1.  and 
monitored  the  conserved  quantity.  Equation  (4.62).  along  the  same  ray.  The  outcome  of 
both  trials  was  qualitatively  similar:  the  computed  conserved  quantity  oscillated  with 
a  period  equal  to  the  interaction  length.  The  difference  between  the  peak  value  and 
the  minimum  value  increased  as  the  grid  size  was  made  larger.  In  addition,  dissipation, 
i.e.  the  drop  of  the  peak  value  as  a  function  of  position  1  increased  as  the  grid  size 
was  made  larger,  and  as  a  result,  the  resulting  local  interaction  length  grew  since  the 
amplitude  of  the  modes  were  attenuated.  While  we  were  unable  to  find  the  cause  for  such 
an  outcome,  we  do  know  that  it  is  not  related  to  the  discretization  of  the  dyy  operator  or 
to  boundary  effects,  since  the  problem  also  arises  in  the  two-dimensional  trial,  which  has 
no  y  dependence.  We  also  tried  changing  the  iteration  discrepancy  tolerance  and  saw  no 
correlation  between  the  value  of  this  parameter  and  the  dissipation.  We  did  find,  however, 
that  the  dissipation  and  oscillation  of  the  conserved  quantities  can  be  made  negligible  by 


130 


Grid  Size  A 

Fluctuation 

4.00 

0.1002 

2.00 

0.0627 

1.00 

0.0168 

0.50 

0.0050 

0.25 

0.0014 

Table  5.1:  Energy  fluctuation  vs.  grid  size.  Equilateral  grid  case. 

making  the  grid  size  small.  We  also  found  that  the  effect  is  much  more  pronounced  when 
A2  =  0  exactly,  which  yields  solutions  with  very  sharp  minimas  in  the  field  variables. 
Table  5.1  shows  the  difference  between  succesive  maxima  and  minima  for  the  second  trial 
as  a  function  of  grid  size,  with  Ax  =  Ay.  We  also  report  the  outcome  of  fixing  Ax  =  0.25 
and  varying  Ay,  in  Table  5.2,  and  the  opposite  settings  are  illustrated  in  Table  5.3.  The 
two-dimensional  trials  for  Ax  =  0.25  and  Ay  =  4  showed  significant  discrepancies  when 
compared  to  the  Runge-Kutta  calculation,  and  the  energy  for  this  case  oscillated  in  a 
somewhat  irregular  pattern.  While  it  is  expected  that  any  discretization  of  the  surface 
system  will  have  inherent  dissipation  due  to  truncation,  especially  manifesting  itself  for 
large  grid  sizes,  it  is  not  at  all  obvious  at  this  stage  of  the  research  that  the  root  cause 
is  truncation  rather  than  some  other  cause. 

To  conclude  this  section,  we  report  the  wall-clock  times  for  three  runs  of  the  surface 
equations,  as  discretized  using  the  Fixed  Point  Method.  The  code  was  written  in  Fortran 
77-because  of  issues  related  to  code  portabilitv-in  a  straight-forward  manner,  except 
that  recursion  was  used  in  the  iteration  procedure.  For  the  size  of  these  runs,  the  use  of 
recursion  was  probably  marginally  slower  than  opting  for  repeated  subroutine  calls.  No 


Grid  Size  Ay 

Fluctuation 

4.00 

0.0018 

2.00 

0.0013 

1.00 

0.0013 

0.50 

0.0012 

0.25 

0.0014 

Table  5.2:  Energy  fluctuation  vs.  Ay.  Ax  =  0.25  fixed. 


Grid  Size  Ax 

Fluctuation 

4.00 

0.1415 

2.00 

0.0028 

1.00 

0.0198 

0.50 

0.0049 

0.25 

0.0014 

Table  5.3:  Energy  fluctuation  vs.  Ax.  Ay  =  0.25  fixed. 
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Machine 

A  =  0.5.  (100  x  100) 

A  =  0.25,  (200  x  200) 

Sun  Sparc  SLC 

7.43 

25.42 

78.8 

Sun  Sparc  2 

2.29 

7.81 

23.13 

Ardent  Titan  2X  Pi 

3.9 

13.9 

44.81 

Table  5.4:  Wall-clock  times  in  seconds  vs.  grid  size  (number  of  grid  points  per  domain) 
for  the  computation  of  the  surface  system  over  the  whole  domain  using  the  Fixed  Point 
Method. 


Machine 

A  =  0.5,  (100  x  100) 

A  =  0.25.  (200  x  200) 

Sun  Sparc  SLC 

0.16 

0.25 

0.50 

Sun  Sparc  2 

0.06 

0.08 

0.15 

Ardent  Titan  2X  PI 

0.08 

0.13 

0.29 

Table  5.5:  Wall-clock  times  in  seconds  for  the  computation  of  the  surface  system  for  all 
values  of  y  at  a  particular  x  using  the  Fixed  Point  Method. 


machine  optimization,  or  floating  point  accelerators  were  utilized.5  The  time  trials  were 
carried  out  with  an  initial  bottom  configuration  /  =  O.Olx.  All  other  parameters  and 
physical  quantities  were  the  same  as  those  used  previously.  The  domain  was  a  square 
with  50  units  to  its  side.  Two  times  are  reported,  the  first  one,  in  Table  5.4,  corresponds 
to  the  total  time  required  to  find  the  field  variables  everywhere  in  the  domain,  and 
a  second  one,  given  in  Table  5.5.  is  the  time  required  to  compute  all  values  in  the  y 
direction,  for  a  particular  x. 


'The  Titan's  vectorizability  was  not  exploited  either.  Otherwise,  its  reported  performance  would  not 
compare  so  unfavorably. 


Chapter  6 

Qualitative  Features  of  the 
Solutions  to  the  Full  Model 

The  main  qualitative  features  of  the  full  model  are  presented  in  this  chapter,  using 
examples  computed  numerically  with  the  Fixed  Point  Method.  The  main  points  of  the 
chapter  are:  To  present  the  effects  of  different  initial  bottom  configurations  and  boundary 
conditions  on  the  surface  and  on  the  eventual  bottom  topography  after  the  passage  of 
many  surface  waves;  and  to  show  that  when  the  slopes  of  the  ocean  bottom  are  very 
mild  and  the  back-wash  negligible,  the  reflected  wave  plays  a  relatively  minor  role  in 
determining  the  shape  of  the  ocean  surface  and  therefore  of  the  sand- ridge  topography. 

6.1  General  Behavior  of  the  Solutions 

To  better  discern  the  effects  of  different  bottom  topographies  on  the  surface  waves  and 
on  the  eventual  bottom  topography  after  the  passage  of  many  waves,  attention  will  now 
be  given  to  the  case  in  which  the  initial  bottom  configurations  are  strictly  z-dependent 
and  the  boundary  conditions  are  constant.  Briefly,  in  this  case,  a  larger  number  of  bars 
form  when  the  gradient  is  slight,  the  distance  separating  the  bars  increases  seaward  for 
the  positively  sloped  case,  and  initial  bottom  discontinuities  in  the  z  direction  tend  get 
“smoothed  out”  after  the  passage  of  many  waves. 
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The  modes  for  waves  that  are  traveling  normal  to  the  shore  over  topography  described 
as  /  =  0.006x  are  displayed  in  Figure  6.1.  Figure  6.2  shows  the  eventual  topography 
of  a  bottom  which  was  initially  the  sloped  but  featureless  profile  of  the  last  example. 
Superimposed,  but  not  drawn  to  scale,  is  the  actual  ocean  surface,  composed  exclusively 
of  an  incident  wave  field  pictured  at  7  =  0.  Figure  6.3  shows  the  eventual  topography 
of  a  bottom  which  initially  had  a  step  in  its  profile.  Note  the  smoothing  effect  due  to 
the  passage  of  many  waves.  All  of  these  figures  had  a  =  0.1,  c  =  0.2,  J  =  0.08,  =  1.8. 

For  the  same  range  of  parameters.  Figure  6.4  shows  the  effect  on  the  surface  and  on  the 
eventual  bottom,  of  an  initial  topography  that  is  approximately  tuned  to  the  interaction 
length  of  the  surface  waves. 


position  x 


Figure  6.1:  <zi  and  02,  for  f{x,y )  =  0.006x.  ai(x  =  0)  =  0.5,  02(1  =  0)  =  0.01 
A  bottom  which  initially  had  gradients  in  the  y  direction  bends  the  water  waves 


position  x 


Figure  6.2:  Ocean  surface  at  T  =  0,  and  below,  bottom  topography  at  T  =  0  and 
T  =  lOOAT.  Not  drawn  to  scale. 

affecting  the  eventual  bottom  topography  by  producing  a  series  of  bars  with  refractive 
features.  Consider,  for  example,  the  case  in  which  the  initial  bottom  topography  is 
f(x,y)  =  O.Oly  ,  with  all  other  parameters  as  before,  except  u>i  =  1.2.  Figure  6.5  shows 
aa  at  T  =  0  and  Figure  6.6  shows  the  refracting  bottom  at  T  =  400AT.  A  striking  way 
in  which  refraction  takes  place  can  be  seen  in  the  case  for  which  the  boundary  conditions 
at  i  =  0  are  y  dependent.  The  case  for  which  f(x,y)  =  0  at  T  =  0  and  the  boundary 
conditions  are  A\  =  0.5  +  O.Oly  and  A2  =  0.1  +  O.Oly,  corresponding  to  an  incoming 
gravity  wave  that  has  slightly  higher  amplitude  at  one  end  than  at  the  other,  is  shown 
in  Figures  6.7,  6.8  for  (^(T  =  0)  and  f(T  =  400AT),  respectively. 

Interesting  configurations  are  achieved  when  the  above-mentioned  effects  are  com- 
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Figure  6.3:  The  fate  of  the  topography  which  initially  contained  a  step,  shown  at  three 
different  times. 

bined.  Figure  6.9  illustrates  the  outcome,  after  T  =  400AT  ,  on  a  bottom  topography  for 
which  -4i  =  0.5  -  O.Oly,  Ai  =  0.1  -  O.Oly,  and  the  bottom  at  T  =  0  was  f(x,  y )  =  O.Oly. 

The  eventual  fate  of  a  bottom  which  initially  was  smooth  but  sloped  in  the  longshore 
direction  is  illustrated  in  Figure  6.10.  The  boundary  conditions  in  this  example  were 
Ai  -  0.5,  Ai  =0.1,  the  initial  bottom  is  described  in  the  figure.  Of  note  is  the  apparent 
growth,  and  motion  of  the  sand  ridges  in  the  shoreward  direction,  particularly  where  the 
water  column  is  deepest. 
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Figure  6.4:  Effect  of  a  tuned  bottom,  /  =  0.5  sin(0.412x)  at  T  =  0,  on  the  eventual 
topography  and  ocean  surface:  Light  solid  line.  Bottom  at  T  =  100AT :  Dark  solid  line. 

6.2  Contribution  of  the  Reflected  Component  to  the  Sur¬ 
face  Waves 

Shown  in  Figure  6.11  is  the  cross-section  of  mode  a\(x,y ),  and  in  Figure(6.12)  a  com¬ 
parison  of  the  eventual  bottom  with  and  without  contributions  from  the  reflected  field. 
Both  figures  were  computed  using  Equation  (4.1),  with  A\  =  0.5,  A 2  =  0.01,  B\  =  0.2, 
and  02  =  0.;  £  =  0.2,  a  =  0.1,  0  =  0.08.  The  bottom  was  f(x,y)  =  0.006x  at  T  =  0. 
The  domain  was  200  units  long. 

As  was  discussed  in  chapter  2,  the  reflected  and  incident  fields  are  completely  decou¬ 
pled,  owing  to  the  assumptions  made  concerning  the  bottom  topography.  The  deforma¬ 
tions  on  the  bottom  topography  due  to  the  reflected  component  are  entirely  determined 
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Figure  6.5:  Refraction  on  the  surface  modes  due  to  the  bottom  topography.  Shown  at 
T  =  0. 

by  the  amount  of  energy  in  the  boundary  conditions.  Hence,  it  is  necessary  to  include 
the  reflected  component  when  the  sea-going  wave  backwash  is  not  negligible. 

If  the  spatial  scales  of  variation  in  the  bottom  topography  in  the  shoreward  direction 
are  of  the  same  order  as  those  of  the  surface  waves,  then  scattering  plays  an  important 
role  in  the  energetics  of  these  surface  waves;  hence  the  reflected  component  must  be 
included  even  if  the  backwash  is  negligible.  The  model  for  the  surface  waves,  in  this 
case,  is  given  by  Equation  (2.75). 
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6.3  Periodic  Solutions  to  the  Surface  System 

To  give  an  idea  of  the  rich  structure  of  the  surface  system,  Equ^iou  (4.14)  is  solved 
in  the  following  examples  using  periodic  boundary  conditions  in  y.  The  parameter  <5 
is  considered  independent  of  frequency.  The  following  graphs  were  generated  using  the 
Fixed  Point  Method,  in  which  the  linear  operator  is  discretized  using  the  Douglas  Scheme. 
The  discretization  yields  a  tri-diagonal  matrix  problem,  with  additional  non-zero  constant 
entries  in  the  upper  right-hand  corner  and  the  lower  left-hand  corner.  This  type  of  matrix 
is  known  as  a  ‘'Jacobi  matrix”,  and  shows  up,  for  example,  in  the  solution  to  the  Toda 
lattice  problem  with  periodic  boundary  conditions.  In  the  figures,  two  periods  in  y 
are  plotted  in  tandem,  the  calculation  being  performed  on  only  one  of  the  two  periods. 
The  domain  has  M  =  240  and  N  =  150,  and  the  the  fundamental  frequency  used  was 
u\  =  1.2.  The  parameters  were  a  =  0.1  and  3  =  0.18.  The  solution  to  the  case  with 
boundary  conditions  A\  =  0. 5+0.1  sin(  ^rx y)  and  Ai  —  0  is  illustrated  in  Figure  6.14.  For 
the  same  parameters,  but  with  the  boundary  condition  A\  =  0.1  sin(^xy),  the  outcome 
is  shown  in  Figure  6.15.  Comparison  of  the  last  case  to  the  case  in  which  6  =  0  is  given 
by  Figure  6.13. 

An  interesting  pattern  arises  in  the  evolution  of  a  case  with  quasi- periodic  boundary 
conditions;  A\  =  0.1(sin(  ^>xj/)+sin(  $-xy)]  and  A2  =  0.  The  outcome  is  Figure  6.16,  with 
the  same  parameter  values  as  in  the  previous  figure,  except  that  6^0.  The  outcome 
shown  in  Figure(6.16)  is  obtained  for  the  boundary  conditions  A\  =  0.1(sin(^)  + 
sin(^)  and  A2  =  0  (i.e.,  quasi-periodic),  and  6  ^  0.  The  numerical  solution  of  this  last 
example  suggests  that  solutions  to  the  surface  system  are  stable  and  periodic. 


Figure  6.6:  Refraction  due  to  initial  bottom  configuration.  Bottom  at  T  —  400A7\ 
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Figure  6.7:  Refraction  due  to  boundary  conditions.  02  at  T  =  0. 


Figure  6.8:  Refraction  due  to  boundary  conditions.  Bottom  at  T  -  400AT. 


Figure  6.9:  Refraction  due  to  boundary  conditions  and  initial  bottom  configuration. 
Bottom  at  T  —  400AT. 


Figure  6.L0:  Evolution  of  bottom 


f.  T=0:  grid.  T=200:  grey.  T=400:  dark. 
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x  position 


Figure  6.12:  Effect  of  a  bi-directional  surface  wave  field  on  the  eventual  bottom  configu¬ 
ration.  Initially,  f(x,y.0)  =  0.006x.  The  dark  line  is  the  bottom  resulting  from  a  strictly 
shoreward-directed  wave. 
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Figure  6.13:  u(x,y)  for  boundary  conditions  A\  =  0.1  sin(^7ry),  A2  =  0,  and  detuning 
parameter  6  =  0. 


Figure  6.14:  u(x.y)  for  A}  =  0.5  +  0. 1  sin(  $iry),  A2  =  0.  and  6^0. 
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Figure  6.15:  u(x,y)  for  Ai  =  0.1  sin(  A2  =  0.  and  <5/0. 


Figure  6.16:  Solution  for  quasi-periodic  boundary  conditions:  A\  -  0.1  [sin(  jrxy)  + 
sin(  $7ry)],  and  A2  =  0.  6  ^  0. 


Chapter  7 

Conclusions  and  Future  Research 
Plans 


This  study  detailed  the  construction  and  implementation  of  a  model  for  the  formation 
and  evolution  of  three-dimensional  sedimentary  structures  on  the  continental  shelf,  based 
on  the  energetic  interactions  of  weakly  nonlinear  long  waves  with  the  shelf  s  sedimentary 
topography.  This  chapter  turns  its  attention  to  the  larger  picture,  discussing  the  main 
conjecture  of  the  study,  as  well  as  the  methodological  aspects  pertinent  to  future  research. 

The  main  conjecture  of  this  study  is  that  a  significant,  but  by  no  means  exclusive, 
agent  for  the  formation  and  evolution  of  longshore  sand  ridges  on  regions  of  the  conti¬ 
nental  shelf  that  are  sufficiently  removed  from  the  shoaling  area  is  the  repeated  action 
of  the  second-order  oscillatory  drift  velocity  that  results  from  the  passage  of  weakly  non¬ 
linear  shallow-water  internal  or  surficial  waves.  The  basis  for  the  conjecture  rests  on 
( 1)  the  close  correlation  between  the  inter-bar  spacing  and  the  length  in  which  energetic 
exchanges  among  the  most  powerful  modes  of  the  shallow  water  waves  takes  place;  (2) 
the  close  correlation  between  the  evolutionary  time  scales  for  the  bars  and  the  time  re¬ 
quired  for  highly  coherent  nonlinear  dispersive  wave  trains  to  impart  sufficient  energy 
into  a  boundary  layer  to  significantly  transform  a  sediment- laden  bottom  topography; 
(3)  the  fact  that  longshore  sand  ridges  are  found  in  areas  in  which  no  wave  breaking 
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occurs  and/or  in  which  the  reflected  field  is  absent  or  negligible;  (4)  the  claim  that  sand 
ridges  with  highly  organized  characteristics  may  be  found  in  regions  in  which  coherent 
weakly  nonlinear  dispersive  waves  exist;  and  (5)  that  the  energy  of  these  waves  is  of  the 
correct  magnitude  to  significantly  affect  the  topography  of  a  sediment-laden  bottom. 

At  present,  neither  the  dynamics  of  sedimentation  nor  those  of  water  waves  are  fully 
understood.  The  model  presented  here  represents  the  conjecture  based  on  current  un¬ 
derstanding  of  both  processes.  If  the  conjecture  is  correct,  the  model  will  improve  in 
predictive  power  as  understanding  of  sedimentation  and  wave  dynamics  improve.  How¬ 
ever,  the  more  important  functions  played  by  the  model  are  that  its  development  viels 
clues  to  ways  in  which  the  conjecture  itself  may  be  refined  and  tested,  as  well  as  providing 
conceptual  spin-offs  (such  as  the  modal  surface  system  that  appears  in  this  study)  which 
are  interesting  independent  of  the  sedimentation  problem  at  hand. 

The  model  in  its  inception  was  two-dimensional.  Based  on  encouraging  comparisons 
with  actual  field  data,  the  three-dimensional  version  was  developed  and  made  the  subject 
of  this  study.  Briefly  described,  the  present  model  couples  a  mass  transport  equation, 
which  controls  the  history  of  the  bottom  topography,  to  a  mathematical  equation,  which 
describes  the  evolution  of  the  most  energetic  modes  of  surface  or  internal  weakly  nonlinear 
dispersive  shallow  water  waves  with  weak  span-wise  spatial  dependence.  In  order  to  solve 
the  coupled  system  one  must  rely  on  the  discrepant  time  scales  of  the  bottom  evolution 
and  of  the  water  waves  to  effectively  decouple  their  interaction,  making  a  solution  by 
iteration  possible. 

In  the  near  future  the  modal  representation  of  the  water  waves  will  be  replaced 
with  a  full  Boussinesq  system,  and  the  effects  of  oceanic  currents  will  be  included  in  the 
model.  Bona  and  Saut  [69]  are  presently  studying  the  different  versions  of  the  Boussinesq 
system  in  order  to  determine,  among  other  things,  which  variant  best  models  oceanic 


waves,  and  which  is  well-posed  as  a  boundary  value  problem.  Additionally,  a  number 
of  issues  brought  up  in  this  study  need  to  be  pursued  to  completion.  These  include 
the  search  for  stable  bottom  configurations  as  predicted  by  the  model,  completion  of 
the  weU-posedness  theory  and  the  Hamiltonian  structure  for  the  surface  system,  and 
development  of  a  stability  result  for  the  iterative  procedure  that  was  used  to  solve  the 
coupled  surface/mass  transport  equations. 

The  sensible  way  to  test  the  conjecture  and  the  model  is,  of  course,  to  examine 
oceanic  field  data.  Comparisons  with  oceanic  field  data  can  assess  the  predictive  pow¬ 
ers  of  the  model:  laboratory  experiments  cannot,  however,  as  they  do  not  scale  well. 
The  task  of  making  field  observations,  particularly  in  the  three-dimensional  case,  is  a 
tedious,  expensive,  and  sometimes  difficult  enterprise.  While  researchers  at  the  INRS 
at  the  University  of  Quebec,  headed  by  Prof.  Boczar  Karakiewicz,  were  able  to  make 
some  comparisons  between  the  two-dimensional  version  of  the  model  and  sand  ridge  data 
[2],  finding  that  the  model’s  predictions  agreed  qualitatively  with  the  height,  spacing, 
and  evolution  trends  of  the  actual  bars,  they  have  not  yet  taken  on  the  task  of  making 
comparisons  in  the  three-dimensional  case.  As  of  this  writing,  the  Quebec  team  is  re¬ 
ducing  field  data  from  the  continental  shelf,  gathered  from  the  ocean  floor  neighboring 
Newfoundland  and  Eastern  Australia. 

There  are,  however,  several  aspects  of  the  conjecture  which  can  be  tested  in  the 
laboratory  as  well  as  in  the  field.  The  drift  velocity  created  by  shallow  water  waves  of 
the  type  identified  here  as  responsible  for  the  formation  of  longshore  sand  ridges  must 
be  observed  and  studied  in  a  laboratory  setting.  Comparisons  between  the  laboratory 
experiments  and  the  drift  velocity  measurements  in  sand  ridge  fields  could  prove  fruitful. 
Additionally,  it  should  be  possible  to  correlate  the  drift  velocity  to  the  shadow  water 
waves  in  question  both  in  the  laboratory  and  in  the  ocean. 


Field  observations  are  needed  to  ( 1 )  determine  the  importance  of  both  the  reflected 
wave  field  and  oceanic  currents  in  determining  the  nature  of  the  drift  velocity  in  sand 
ridge  areas;  (2)  correlate  in  some  way  the  beginning  and  end  of  ridge  fields  and  the 
physical  location  at  which  water  waves  are  created  and  eventually  destroyed;  (3)  track 
the  relevant  wave  spectra  in  order  to  see  evidence  of  the  predicted  pattern  in  energetic 
interaction  lengths  and  its  correlation  to  features  of  the  bottom  topography;  and  (4) 
determine  what  sort  of  sand  ridge  configurations  are  stable  and/or  non-migratory. 

Laboratory  observations  are  required  to  determine  how  well  the  various  Boussinesq 
systems  model  the  weakly  nonlinear  shallow  water  waves  and  to  confirm  the  existence  of 
recurrence-like  solutions  over  long  propagation  lengths.  Additionally,  more  experiments 
aimed  at  furthering  our  understanding  of  the  motion  of  sediment  in  the  boundary  layer 
are  needed. 

Computational  experiments  are  currently  being  planned,  aimed  at  exploring  the  na 
ture  of  recurrence-like  solutions  in  nonlinear  dispersive  equations,  such  as  the  Boussinesq 
equation;  other  experiments  will  explore  the  stability  and  interdependence  of  the  trun¬ 
cated  modal  solutions  to  these  equations. 

In  conclusion,  this  study  has  produced  a  wealth  of  interesting  and  fundamental  ques¬ 
tions.  While  comparisons  between  field  data  and  the  two-dimensional  model  are  very 
encouraging  and  this  three-dimensional  extention  should  therefore  find  applicability  in 
the  real  world  environment,  any  topographical  chart  of  the  continental  shelf  provides  a 
good  reminder  of  the  long  path  yet  to  travel  toward  a  complete  understanding  and  model 
of  the  full  problem.  If  this  study  has  piqued  the  curiosity  and  compelled  the  reader  to 
take  a  closer  look  at  sandbars,  it  will  have  succeeded. 


Appendix  A 


Higher  Order  Theory 


Considered  here  is  the  second  order  contributions  to  the  proposed  model.  This  section 
is  included  as  an  appendix  because  of  it’s  exploratory  nature.  The  value  of  these  cal¬ 
culations  resides  in  the  possibility  of  discerning  if  any  fundamentally  new  contributions 
may  arise  from  a  careful  inclusion  of  these  higher  order  terms.  The  very  tedious  process 
of  generating  the  surface  contributions  at  this  order,  and  the  daunting  problem  of  for¬ 
mulating  the  drift  velocity,  is  a  veritable  project,  even  for  a  symbolic  solver.  The  higher 
order  expressions  were  derived  as  ca'efully  as  possible,  nevertheless,  it  is  very  possible 
that  algebraic  errors  were  made. 

The  irrotational  condition  to  this  order  is 


a2  :  u2y  -  v2l  -  vix  =  0. 


(A.ll 


and  the  continuity  condition. 


0,2  ;  l2t  +  a2x - -j—  =  F2 {rjo.  U0.  vq,  t)\,  ui,  t’i.  G\  x,  X,  y,  t), 


(A.2) 


where  G(X.  y,  T)  =  ^V'r)  t  an(j 
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E2  =  ~V\y  + 


a2 

y  ^lyyt 


—  WiX  —  Ml  T)qx  -  GuXl  -  TftUox  + 


20  {nixXt+^Vl  lit  ) 


—  uoVix  —  7ouix  —  Gyi'o  —  G  \  uo  —  i'oHoy  -  Gi'oy  —  rjQi'Qy  +  —  L.JVa0ns 

0  VOXXt  40  GxVO-ct  ,  40  Gr) QiX_t 


(A. 3) 


-UoVoX  -  Guqx  -  noUox  H - 3 

Reiterating  Equation  (2.60). 


a2:  Ct]2  =  ^2(??o.wo,Wo-'?i,  «i,  0\,G\x,  X,y,t)  (A.4) 

where  the  linear  operator  is 

C  =  dtt-dIX-  2!3£££i,  (A. 5) 

and  the  inhomogeneous  term  is 

C>2  =  (1  +  02dtt/3)f]iyy  +  G{  1  +  '202dtt/3)T)iXI  +  2(  L  +  i32dtt/3)V\xX 
( 1  +  02dtt/3)T]oxx  +  G(  1  +  2i32dtt/Z)V0yy  +  2(  1  +  02dtt/3),noxx 
+Gxi  1  +  ■id2dtt/3)  +  rfox  +  GyTioy  +  '202GyyTiott/3  +  402GyT)otyy/3  (A. 6) 

~(m  «0  +  7o«i)r(  +  (uiUo)xx  +  (u2o)xX  ~  ( u0Tjo)xt  +  G{ul/2)xx 
4-(u0/2)yj/  +  (  Vq/2)xx  +  (  0o/2)yy  ~  (  r)Otvo)y- 
Applying  the  compatibility  condition,  we  get,  after  a  substantial  amount  of  algebra. 


Air  +  ie f D\E\A\  —  iaF\A\yy  —  iotD\S\t  ‘^x{A\a.2  +  aj A2) 
A2X  +  iefD2E2A2  -  iaF2A2yy  +  i2aD2S2e+^xa\ Ax 


ni(x.j/,a1,a2) 

ft2{x,y.ax,a2). 

( A. 7) 


for  the  equations  of  the  modes.  The  constants  on  the  left  hand  side  are  given  by  Equation 
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(2.73),  and  the  inhomogeneous  terms  are 

=  -piaixyy  -  p2f(x,y)aiI  -  p3alxT  -  p4fA^-y)ai  +  p5fyy(x,y)ai 

+p6fy(x,  y)aiy  +  p7/(x,  y)aXyy  +  {-p8a2aJJ  +  pga2xa\  +  pxof(x,y)a\a2 
+pn(aja2)x  -  pu(a^a2)yy  -  pX3(a\ya2)y  -  pu(a\a2y)y  -  pna2ya\y}e~‘dx 
+Pi6(aIai  +  2axa2a2) 

(A.8) 

and 

^2  =  -pna-2xyy  -  Pi8f{x,y)a2x  -  piga2xx  -  p2xfx{.c, y)a2  +  p22fyy(x,y)a2 

+P23fy(x^y)a2y  +  P24f{x,y)a2yy  +  { -p2SQl Ojx  +  P2ef{x,  y )<*]  +  p27( of  ) 

P2s(  )yy )  +  P29(aiy)2  +  P3o(aiya1)y}e+l6z 

+P3i(a2a2  +  2axa^a2). 

(A. 9) 

The  coefficients  appearing  in  the  fi  terms  above  are  explicitly  given  below: 


Pi  =  2  F2 

p2  =  2eD\F\E\l  a 

P3  =  ifi/a  (A.  10) 

4/?2u;? 

P4  =  £2^(1--^) 

2  2 

Ps  =  uxDi 

p6  =  4isFx  D\k\  + 

p~  =  4ieFfDxEi 

ps  =  4DiF2(k2  -  ki)(u22  -  u>i) 

P9  =  4D\Flu>2(k2  —  k\)(u)2  —  uji)/sj\ 

Pio  =  -  ki  )1[2DiEx(k2  -  k\ )  -  kx] 
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Pn 
Pi  2 
Pl3 
Pi  4 
Pi  5 
Pl6 
Pi  7 
Pl8 
Pl9 
P21 
P22 
P23 
P24 
P25 
P26 

P27 

P28 

P29 

P30 


{4Fi(*2  —  ^i)[2£iFi(^2  —  &1  )  —  ^'l]  —  (<*2  —  ^1  j(“J2/^2  +  “f'l/fcl  )} 

u2i 

AiaFxF2Dx^k2(k2  -  kx)  ,  ,  , 

"l - Tn - *1*2  ”  ^1^2] 


^1 


2D, 


AiaF\F2D\uj\lui\ 
AicxF,  F2D,uJ\ 


AiaF\F2D\(k2  -  k\  )2 /ui\ 
otD, 

2F\ 

2zD2F2E2j  a. 


3iF2/au2 

z2D2F2(  1  — 
2  , 

l£~PiijJ2D2 


Ai32ul 


)/a 


AUF*  D2k2  +  "zS0^>jJ2D2I  k2 
o 

4  ieF,  DXEX 
i8cJi  D2F\ui2 
2 kD2[^-  -  1  }kl/uj2 

2D2[-p--2u\Fx-k,] 


D 


~^4m(i+4)] 

u>2  D  2  2  1 

2iaD2/u>2 

iaJ\D2/u2k\  )2 

a£)2. 


P31 


(A. 11) 
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The  drift  velocity  to  second  order  in  a  is 


<U2)  =  (u2)  +  (J  u,di-Vuo>  +  { J  j  uodt'VuodtVuo) 

/t  1  ft  -  ft 

uodi-Vui)  +  -({J  u 0di}T-Hnn- u 0dt). 


(A.  12) 


which,  after  weak  y  dependence  scaling  is  adopted,  can  be  expressed  in  component  form 


U  —  (u2)  +  (fl  vodtuoy)  +  { ff  uidtuQx)  -+■  (f‘  w\diuon)  +  (f*  uodtuix)  +  ( J‘  wodiu\n) 
{{f  f  u0dt'u0xdi  +  f‘  fl  wodt'u0ndi}uox+ 

{f‘  Jf  u0dt'w0Tdi  +  f‘  f i  w0dt'w0ndi}u0n+ 
i(/'  U0di)2U0xx  +  j(f‘  IL'odtfuonn) 

V  =  (v2)  +  ( f‘  vodtv0y)  +  (f‘  uidtvox)  +  (J*  Wjdtvon)  +  { f ‘  u0div\x)  +  (/'  u'odivln) 

({fl  ft  uodt'uoxdt  +  f‘  wodt,Uondi}i\ix+ 

if‘  fl  Uodt’woxdt  +  /*  ft  W0dt'tl'0ndi}V0n  + 

]2{  fl  U0dt)2V Oxx  +  W0di)2V0nn). 

(A. 13) 

The  longshore  drift  velocity  was  calculated  for  which  the  second  order  velocity  in  the 
boundary  layer.  Equation  (3.21),  is  needed.  The  calculation  was  carried  out  by  isolating 
the  contributions  at  the  first  and  second  harmonics  to  the  velocity.  Once  the  velocity  was 
calculated,  Equation  (A.  13)  was  computed  explicitly.  Finally,  U  gets  integrated  over  the 
depth  of  the  boundary  layer  since  the  density  distribution  for  the  sediment  is  assumed 
constant.  The  resulting  expression  for  the  longshore  mass  transport  velocity  is 


158 


(A.  14) 


where  M,  X,  P ,  L  are  complicated  coefficients  that  depend  on  the  frequency  and  wavenum¬ 
ber  of  the  waves,  the  boundary  layer  thickness,  and  the  parameter  J.  Finally,  we  add 
the  above  contribution  to  the  mass  transport  equation,  which  now  looks  like 


d-~~—  =  —[Mix  +  a(^2x  +  )]•  (A. 15) 

oT  pQ 

An  illustration  of  the  contribution  of  the  second  order  theory  is  shown  below.  Figure 
A.l  shows  the  cross  section  of  the  higher  order  contributions  to  the  surface  wave  field, 
and  Figure  A. 2  gives  a  comparison  of  the  eventual  bottom  topography  which  is  formed 
by  the  action  of  the  higher  order  drift  velocity  to  the  case  in  which  no  higher  order 
contributions  are  used.  The  bottom  initially  had  a  profile  f(x)  =  0.006,  the  frequency 
was  u>i  =  1.8.  The  parameters  were  set  at  a  =  0.1,  ^  =  0.08,  s  =  0.2. 

Figure  A.l  and  Figure  A.2  display  the  higher  order  effects  for  the  same  problem  and 
parameters  used  to  generate  Figure  6.1  and  Figure  6.2.  Although  the  plots  for  the  higher 
order  contributions  have  been  produced  with  an  exagerated  vertical  scale,  it  is  not  clear 
whether  in  reality  ,  these  quantities  do  indeed  eventually  diverge. 
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Figure  A.l:  Higher  order  contributions  to  the  surface  wave  field  at  T  =  0.  The  vertical 
scale  has  been  exagerated. The  lower  curve  represents  the  bottom.  A\: - ,A2:- 
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Figure  A.2:  Higher  contributions  to  the  bottom  topography.  The  bottom 
f(x)  =  0.006.  Eventual  bottom  with - ,  and  with  no . 


was  initially 
higher  order 


contributions. 


Appendix  B 


Slightly  Interacting  Resonant 
Quartets 


This  appendix  contains  the  expressions  for  the  lowest  order  surface  wave  modal  ampli¬ 
tudes  for  the  case  of  quartet  interactions. 

The  presentation  is  limited  to  the  incident  wave  field.  The  relation  among  the  fre¬ 
quency  and  wavenumbers,  that  Uj  =  juq,  k2  =  2k\  -  6.  and  k3  =  3ki  -  A  is  given  by 
the  dispersion  relation.  The  procedure  is  the  same  as  the  two  mode  case.  Substituting 
Equations! 2.65)  and  (2.67),  with  j  =  1.2,3,  into  the  compatibility  condition.  Equation 
(2.71)  yields  the  following  system: 

a\i  +  itfD\E\a\  -  ia.F\a\yy  +  iQD\S2iie~'^za\a2  +  iaD\S32iel^xa2a3  =  0 
a2x  +  iefD2E2a2  -  iaF2a2yy  +  iaD2S2e+lf>za\  +  iaD2Szi2e'^z a*^  =  0 

a3x  +  iefDzEzaz  -  iaF3azyy  +  iaD3S3e~l^za\a2  =  0, 

(B.l) 
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to  0(6/ X).  The  constants  are 


32  u>* 

D}  =  1/2(1 - 3^ 

2 

Ej  = 


=  1/2  kj 

s  =  k±kj 


53  =  *^{*2 +  *1+3^1^  +  ^)} 

52  = 

■S*.,/  =  '1u^{kx  -  kj  +  I^J i(xf  +  ^-)} 

The  boundary  conditions  are  similar  to  Equation  (4.2).  except  that  there  are  three  modes 
rather  than  two  that  need  to  be  pinned  down  at  the  boundary.  Thus 


a.j(0,y)  =  Aj, 


with  j  =  1,2,3.  plus  appropriate  boundary  conditions  on  the  lateral  sides  of  the  domain. 
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